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Abstract 


Identifying the common 3-D substructure in two drug or protein molecules is a 
important problem in synthetic drug design and molecular biology. This probler 
can be represented by the following geometric pattern matching problem : give 
two point sets A and B in three-dimensions, and a real number e > 0, find tb 
maximum cardinality subset S C A and an isometry Z, such that each point < 
1(S) is within e distance of a distinct point of B. Since it is difficult to soh 
this problem exactly, in this thesis we have theoretically studied this problem an 
proposed several approximation algorithms with guaranteed approximation ratio. 

Our algorithms can be classified into two groups. In the first we extend tt 
concept of approximate decision algorithms for approximate congruence of poii 
sets in 2-D, to approximately maximize the size of S. All the algorithms in th 
class exactly satisfy the constraint imposed by e. In the second class of algorithn 
this constraint is satisfied only approximately. In this case, we improve the existir 
approximation ratio for this class of algorithms, while keeping the time complexr 
unchanged. For the existing approximation ratio, we also propose algorithms wii 
substantially better running time. We also suggest several improvements of the bas 
algorithms, all of which have a running time of 0(n 8 ' 5 ). These improvements main 
consist of using randomization, and exploiting some general structural properties 
protein molecules. As a result, our final algorithm has a running time of 0(n^ logn 
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Chapter 1 
Introduction 


1.1 Background 

Determining the structural similarities between a set of molecules is a central issue 
in both synthetic drug design and in studying biomoleculax recognition and interac- 
tion of proteins. The need to solve this problem arises because of two assumptions 
in computational chemistry - (i) drug activity is obtained through the molecular 
recognition and binding of one molecule, called the ligand, to the pocket of another 
molecule, called the receptor [Boy90] (ii) functional properties of proteins depend 
on some typical parts of their three dimensional structure, called 3-D structural 
or binding motifs [BT91]. Thus two drug molecules having a large common sub- 
structure might show similar drug activity because of the substructure binding or 
docking itself into the pocket of the receptor. On the other hand, since families of 
proteins retain a common underlying 3-D structure, identification of this is crucial 
in understanding the mechanism by which they work. The receptor-ligand recogni- 
tion and binding is also encountered in a large number of other biological processes, 
such as cell-cell recognition, enzyme catalysis (and inhibition), and in regulation of 
development and growth. 

Again, in the case of structure-based drug design, since the geometric structures 
of a very few receptor molecules are known, trying to develop pharmaceutical drugs 
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for these receptors is an important research problem [BGSS, MBD + 93]. Finding 
out the common substructure that is contained in some active conformation of each 
of the ligands, that have been experimentally found to interact with the considered 
receptor, can give important insight into the structure of the receptor. Since the 
common substructure probably docks into a pocket of the receptor, the shape of the 
pocket is the geometric complement of the substructure, a determination of which 
can lead to the design of more effective pharmaceutical drugs. 

Other applications of this problem are in (i) molecular database screening, where 
we might want to detect a structural pattern in a large set of structures, or retrieve 
functionally equivalent, but structurally novel molecules from the database [Wil95], 
and (ii) in suggesting alignments of molecules for input to CoMFA (Comparative 
Molecular Field Analysis ) and other 3-D QSAR (Quantitative Structure-Activity 
Relationship) methods [BGSS, ITY + ]. 


1.2 The Common Substructure Identification Prob- 
lem 

Considering the applications detailed in the last section, we address the following 
problem : given two molecules, find the maximal common rigid sub-unit contained 
in both the molecules. This is equivalent to determining the rotation and tr anslat ion 
of one molecule, relative to the other, such that there is a large fit or superimposition 
between the two molecules. Once such a common substructure is isolated, it can 
be used to detect its presence in other molecules exibiting similar activity, and can 

also be chemically evaluated to identify the degree to which it is responsible for the 
activity in question. 

In the receptor-ligand binding problem, the receptor molecules mostly being 
proteins are fairly large. The ligands when proteins, are of the order of several 
thousand atoms, and when chemicals or drug molecules, are in the range of ten 
to a few hundred atoms. Thus to systematically deal with such large volumes of 
spatial data, the need for efficient algorithms comes into picture. Although there 
is a considerable volume of literature which address this problem, they still do not 
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seem to be sufficient from a viewpoint of computation time [Aku96], and hence leave 
scope for further research. More importantly, almost none of the known methods 
provide a theoretical guarantee for the quality of the obtained results [Aku96] . So we 
study this problem from a theoretical point of view and propose algorit hms which 
have guaranteed performance bounds in terms of the quality of the results. Towards 
this end, we view this maximal common substructure identification problem in an 
abstract geometric setting where each atom is considered to be a point in 3-D space. 
Thus a molecule is treated as set of labeled points in 5ft 3 . This reduces to finding 
the point set of maximal cardinality that is contained in both the given point sets, 
representing the two molecules. Now, since atom positions are fuzzy, it is impractical 
to consider an exact match between two points. Hence considerable tolerance is 
allowed and two points pi and p% are said to match whenever \pi — p 2 1 < e, where e 
is a predefined constant, usually referred to as point location error. 


A protein may be viewed as a sequence (linear chain) of amino acids which folds 
to generate a compact domain with a specific three-dimensional structure. Many 
substructure matching algorithms for proteins, making use of this sequence property, 
are based on character string comparison algorithms. However such algorithms can 
hardly find useful substructures, whose nature is intrinsically 3-D [PA94, FNW92]. 
So we make no assumption about the linear ordering of amino acids in a protein 
molecule. Again, although from a chemistry point of view, the problem with drug 
molecules might be different from that of proteins, from a geometric standpoint both 
the problems are similar, except for the fact that protein molecules are usually much 
larger than drug molecules. Since in this thesis we will view the problem only from 
a geometric point of view, henceforth we will not distinguish between the two cases 
and will be concerned only with the abstract problem, except in Chapter 4 where 
we suggest some improvements in our algorithms, making use of properties specific 
to protein molecules.. 
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1.3 Geometric Pattern Matching 

Geometric pattern matching is a well studied problem in computational geometry, 
having applications in computer vision [HH92] and computer graphics. In the sim- 
plest case, given two finite point sets P and Q , it is required to be found if they are 
congruent, i.e. if there exists a transformation T consisting of translation, rotation, 
and possibly reflection, such that T(P) = Q. O(nlogn) algorithms to solve this 
problem in 3-D were given by Atkinson [Atk87] and also by Akutsu [Aku92j. For 
any d > 3 and point sets in 5R d , an 0(n d ~ 2 logn) algorithm was given by Alt et al. 
[AMWW88]. A generalization of this, called the congruent copy detection or CCD, 
asks if point set Q is congruent to some subset of P [CGH + 93, GM094, dRL95]. 
However, the abstract problem that we formulated in the last section is more sim- 
ilar to a further generalization of the CCD, called the LCP or the largest common 
point set problem [ATT97, AH]. In this we ask for a set R of maximum cardinality 
that is congruent to some subset of P and to some subset of Q at the same time. 
Such an R is called the largest common point set between P and Q. With |Pj = n, 
|Q| = m and n > m, an 0((n 1A3 m L77 + n 2 ) logn) algorithm for LCP in 2-D and an 
0((min{n L8 m 3 + n 3 , n 4,7 }) logn) algorithm for LCP in 3-D were given in [ATT97]. 
The problem of finding LCP with congruence replaced by similarity was considered 
by Irani and Raghavan in [IR96]. 

In all the above cases, input data is assumed to be exact, free from any noise 
or inaccuracy; and thus only exact matching between two points was considered. 
However, since such an assumption is not realistic in practice, the problem of ap- 
proximate congruence or e-congruence was formulated by Alt et al [AMWW88]. 
Under this formulation two points are said to match whenever they are within e 
distance of each other. All the problems mentioned in the last paragraph have their 
corresponding analogs under e-congruence. Given two 2-D point sets, an algorithm 
to decide wheather they are e-congruent under general isometry with Eucledian 
metric, was given in [AMWW88], with a running time of 0(n 8 ) (in contrast to the 
O(nlogn) algorithms [AMWW88, Ata84] for testing exact congruence). This deci- 
sion problem is easier if the correspondence between the matching points is known in 
advance. Alt et al. [AMWW88] and Iwanowski [Iwa90] gave O(nlogn) algorithms 
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for some restricted decision problems and Imai et al. [ISI89] gave 0(n z logn) algo- 
rithms for the general case with Eucledian metric and 0{ri) algorithm for maximum 
metric. Similarly it was shown by Arkin et al. [AKM + 92] that improved running 
times can be obtained if the e-balls around each point are disjoint. 

A related optimization problem is that of finding the minimum tolerance value 
e, such that the two point sets are e-congruent [AMWW88, GM094]. A similar 
problem which has been extensively studied is that of finding the isometric transfor- 
mation that minimize either the directed or undirected Hausdorff distance between 
two point sets [CGH+93, CK92, HK90, HKK92, HKS91, Rot91]. 

However in our case, which is thus finding the LCP of two point sets under e- 
congruence, the value of e is fixed and a point of one set is mapped to only one point 
of the other set, i.e. we require a bijection. This is not necessarily true in the case 
of measuring resemblance by the Hausdorff distance, where the mapping is defined 
by associating a point with its closest neighbor in the other set. 

1.4 Approximation Algorithms for Approximate 
Congruence 

Almost all algorithms for approximate congruence, except for the very trivial cases, 
suffer from high running times, even in two-dimensions. Indeed, Rucklidge [Ruc93] 
gives evidence that such methods must have high running times. This motivates the 
need for approximation algorithms for approximate congruence, which might have 
better time complexities. Note that the term approximate is used twice with two 
different senses, the first means that the algorithms are not exact in terms of the 
solutions they output, and the second refers to a property of the point sets. 

1.4.1 Two Approaches for Designing Algorithms for Ap- 
proximate Congruence 

There are two main approaches towards this direction. The first is to specifically 
design an approximation algorithm, taking approximate congruence into considera- 
tion. The second is to take an algorithm for the exact matching model and use its 
approximate version, where approximate equality is used instead of exact equality, to 
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test matches. Although there is no marked distinction between the two approaches, 
in the available literature mostly the second approach is taken and is also probably 
preferred in practice. However as already mentioned, it is difficult to bound the 
performance of such algorithms. Work that falls into this category include that due 
to Irani and Raghavan [IR96] and also the recent paper by Finn et al. [FKL+97]. 
Exceptions to this approach is that of Goodrich et al. [GM094] and also the work 
due to Schirra [Sch92, Sch88], Heffernan [Hef91] and Akutsu [Aku96]. 

1.4.2 Algorithms with Guaranteed Performance Bounds 

As mentioned in Section 1.2, we theoretically study the common substructure iden- 
tification problem, which with our present formulation is to find the LCP of two 3-D 
point sets under e-congruence, where e is the given point location error. Towards 
this goal we propose approximation algorithms which fall into the first of the two 
categories mentioned in the last subsection. The methods which we use for designing 
these algorithms are based largely on existing techniques (mainly for 2-D point sets 
of equal cardinality, and for testing the congruity of the entire point set), with a 
number of adaptations for our particular problem. Our problem is a generalization 
of the problem of testing the e-congruence of . two point sets in 2-D, both in terms of 
the dimension of the point sets and in terms of the conditions for congruity. Hence 
it is atleast as hard as that problem, and most likely much harder. An 0(n 8 ) al- 
gorithm for testing the e-congruence of two 2-D point sets was given by Alt et al. 
[AMWW88]. The approximation algorithms that we propose have time complexity 
close to this. 

I (a, / 3) -Approximate Algorithms 

One class of approximation algorithms for testing e-congruence, developed by Schirra 
[Sch92, Sch88], Heffernan [Hef91] and Heffernan and Schirra [HS92] solves the deci- 
sion problem for values of e not too close to €<**, where e^ is the smallest tolerance 
value for which the two point sets are e-congruent. If c is too dose to e^t then 
the algorithm gives up and returns DON’T KNOW , otherwise it solves the problem 
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correctly. These algorithms are called approximate , since they are not always able 
to solve the problem. For equal cardinality point sets, an algorithm is called an 
(a, /3)-approximate algorithm for approximate congruence if it solves the decision 
problem for all values of e outside the interval [e opt — a, € opt + /3). 


I Violating the e-Constraint 

Considering the optimization problem of finding the minimum tolerance value of e 
such that two point sets are e-congruent, Goodrich et al. [GM094] gave an approxi- 
mation algorithm which finds a rigid transformation, which enables 5eop r congruence 
between the two point sets, where e^t is the required minimum value and S > 1. 
For 3-D point sets and general isometry, 5 = 8. 

Although it is not clear how to formulate the LCP problem in this model, since 
now there are two quantities to optimize (the subset size and the distance between 
the individual points) [ATT97] , Akutsu [Aku96] was able to make use of this tech- 
nique to obtain approximation algorithms for the protein structure alignment prob- 
lem. In terms of our problem, i.e. finding the LCP of two 3-D point sets under 
e-congruence, this amounts to satisfying the e-constraint approximately. His algo- 
rithm returns a subset of size larger than the largest common point set that can 
be obtained under e-congruence. However, each point of such a subset is atmost 8e 
distance away from the corresponding point in the other set. 

To remedy this situation, a technique was proposed by which the factor of 8 can 
be reduced to any arbitrary constant 5 > 1 , with only increasing the time complexity 
by a constant factor. However, the resulting algorithm is not practical because the 
constant factor depending on 5 is very large, 

l Our Work 

Our approximation algorithms are based on the last two classes of algorithms we 
have just mentioned, and are motivated by the papers due to Akutsu [ATT97], 
Schirra [Sch92] and Goodrich et al. [GM094]. The factor 8 in Akutsu’s algorithm 
being quite large, may not be tolerable in practical situations. We give an algorithm 
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where this factor is reduced to 2, without incurring any large constant factors. 
Next, instead of approximating the constraint imposed by e, we propose algorithms 
which approximate the size of the largest common point set, and give upper and 
lower bounds on its size. Our algorithms are based on the ideas introduced by 
Schirra [Sch92], which were recently also used by Dehne et al. [DG97] to solve the 
unrestricted point set stereo matching problem. However they are nontrivial gener- 
alizations of [Sch92], because Schirra’s algorithms are for solving the e-congruence 
decision problem for two equal cardinality point sets and are thus able to make use 
of the centroids of the point sets. Our approximation algorithms for finding the LCP 
have subtle but important important differences with [Seh92] because they involve 
an optimization aspect and moreover cannot make use of the centroids of the point 
sets. 

We also show that making use of randomization can reduce the time complexities 
to a large extent. This technique is now fairly standard for this class of problems and 
was recently also used by Finn et al. [FKL + 97] for the common substructure identi- 
fication problem for drug molecules. However as already mentioned in Section 1.4.1, 
their approach is not systematic and does not give any performance guarantee on 
the size of the common subset. Lastly we show that when the molecules whose 
substructure to be identified are protein chains, then by exploiting some general 
structural properties of proteins, the time complexity of our resulting algorithms 
reduce by a considerable extent. 


1.5 Organization of the Thesis 

In the next chapter we briefly outline the important approaches which are popularly 
used for the substructure identification problem and review related work. In Chapter 
3, we formally state our abstract geometric problem and give various approximation 
algorithms. Chapter 4 is concerned with improvements of these algorithms, mak- 
ing use of randomization, the underlying geometric nature of the point sets, and 
the structural properties of protein molecules. Chapter 5 contains conclusions and 
directions for further research. 


Chapter 2 


Basic Approaches 


As we have mentioned in the last chapter, all algorithms for the approximate con- 
gruence problem tend to have high running time. For example, the algorithm to test 
the e-congruence of two point sets in 5ft 2 has a running time of 0(n 8 ) [AMWW88]. 
Further, such algorithms are difficult to implement, as they need to compute intersec- 
tions of complex algebraic surfaces [IMV]. On the other hand, from the viewpoint of 
practical implementation, algorithms for exact congruence with approximate equal- 
ity instead of exact equality to test matches, perform better in terms of ease of 
implementation and running time [IMV, ATT97]. Hence this is the preferred ap- 
proach in practice. 

In this chapter we outline three main approaches popularly used for solving the 
common substructure identification problem. These are alignment [HU90], voting 
[ATT97],. and geometric hashing [Wol90]. Although actually formulated for the 
exact matching model, all of them, specially the last two, for aU practical purpose 
can be used for real world data. In all the cases we mention both deterministic and 
randomized algorithms based on these schemes. Finally we briefly review related 
work in this field, mostly those pertaining to the protein structure alignment problem 
and which have made use of the structural properties of the molecules in some form 
or the other. Ail the algorithms presented here are for point sets in 3-D space. 
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2.1 Alignment 

Given two point sets P and Q, in the deterministic version of the alignment scheme, 
for every triplet of points ( pi,p2,Pz ) from P, all triplets of Q are found 

such that the triangle A P = (pi,P2,P3) is congruent to the triangle A q = (<?i, <?2, 53)- 
Now for each Aq, the transformation 7 pq is calculated so that Tpq(Ap) = A q. This 
transformation Tpq is now applied to the entire set P and the size of the intersection 
Tpq(P ) n <3 is counted. The transformation for which the size of this intersection 
is maximized is the required transformation and gives the largest common point set 
that is contained in both P and Q. Let us denote it by LCP(P,Q). 

There are a number of ways in which we can reduce the running time of this 
algorithm, however, at the cost of a small failure probability, using standard ran- 
dom sampling procedures. Firstly, instead of using the entire point set P, we can 
randomly sample a subset R of P and use all possible triplets of R, instead of P, to 
calculate the transformations. If R contains atleast three points of LCP(P , Q) then 
this procedure successfully finds it, with running time less than that required in the 
previous case. To further reduce the running time, we randomly sample a subset R 
of P and a subset S of Q. For all triplets of points in R and all congruent triplets 
of points in S, the corresponding transformation is calculated. This transformation 
is now applied to the entire point set P, to find out the number of matching points 
in Q. If |LCP(P, Q) n R n T -1 ^)! > 3, where T is the transformation for which 
T(P) D Q is maximized, then this scheme successfully finds LCP{P,Q). For this 
the sizes of R and S are so chosen that \LCP(P, Q) n R n T“ 1 (5')| is atleast three 
with high probability. 

To extend this algorithm to incorporate e-congruence, two points are said to 
match whenever they are within e distance of each other. In this case, for every 
triangle A P = (pi,P2,P3), all triangles Aq = (71,52,53) are considered for which , 
d(Pi,Qi) <€,(* = 1,2,3). For calculating the transformation Tpq, since there can 
be numerous transformations which map a point pi into the e-ball around the point 
5», usually some heuristic is used to resolve the situation. This transformation is as 
usual applied to the entire point set P to calculate the number of matching points 
in Q. The transformation which yields the maximum number of matching points 
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is used to calculate LCP{P,Q). Now, because of the heuristic that was used to 
compute the transformation, it becomes difficult to give a theoretical guarantee of 
the size of the common point set obtained, with respect to the actual result. However 
in a practical situation it suffices to use this scheme. 

* 

2.2 Voting 

For pairs of points (pi,p 2 ) € P x P and (qi,q 2 ) € Q x Q, with d(pi,p 2 ) = d{q r ,q 2 ), 
and any transformation T, let rnult PlP2!qiq2 (T) be the number of pairs x € P and 
y € Q satisfying T( A P ) = A Q where A P and A Q are the triangles (pi,p 2 ,x) and 
(qi,q 2 ,y) respectively. For point sets P and Q as before, if \T{P) nQ| = k, it is 
easy to see that rnult PlP2tqiq2 (T) = k — 2 for any p x and p 2 belonging to T(P) PI Q 
and gi = T- 1 (p 1 ),g 2 = r- 1 (p 2 ). 

If \LCP(P,Q)\ = K, then clearly max plP2)9l92 (max r mult PlP2>qiq2 (T)) = K - 2, 
and the transformation achieving this maximum is the transformation that gives 
LCP(P,Q). For pairs (pi,p 2 ) € P x P and (?i,g 2 ) eQxQ, max r mult PlP2 , qiq2 (T) 
can be computed by considering all the congruent pairs of triangles A P = (pi,p 2 , x) 
and A = (q x , q 2 , y), x € P and y eQ, and letting each such pair cast a vote to the 
■transformation T for which T(A P ) = A q. Let us call this process local voting for 
the pairs (pi,p 2 ) and (q x , q 2 ). In the deterministic version of the algorithm, this local 
voting process is executed for all possible pairs (px^) and (qi,q 2 ), Pi € P,<fc € Q 
(i = 1,2), and the best result is taken. The transformation receiving the maximum 
number of votes is used to compute LCP(P, Q). 

Now since in the above deterministic algorithm, the transformation T receiv- 
ing the largest number (K — 2) of votes, does so for every possible pair (pi.p^, 
(T(pi),T(p 2 )) such that pi,p 2 G LCP(P,Q), there is a lot of redundancy in this 
process. This redundancy is difficult to avoid deretministically unless we know the 
value of K. Randomization can be used to address this problem, but again at the 
cost of a small failure probability, in exactly the same way as was done in the case 
of algorithms based on alignment. Firstly we can randomly sample a subset RoiP 
and do a local voting for all possible pairs (pi , p 2 ) € RxR and (q x , q 2 ) € QxQ. Here 
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once the pairs (pi,p 2 ) and (q u q 2 ) are fixed, the local voting is done exactly as in the 
deterministic case, scanning the entire sets P and Q (not restricting to R). Clearly 
if \Rn LCP(P,Q)\ > 2, this scheme successfully finds out LCP(P,Q), but with a 
running time less than that required in the deterministic case. To further reduce 
the running time, we can as in the case of alignment, randomly sample subsets R 
of P and S of Q and do a local voting for all possible pairs (pi,p2) € R x R and 
(?i, <h) € S x S. LCP(P, Q ) is successfully found if \LCP(P,Q) n JRnT~ l (5)| > 2, 
where T is the optimal transformation. 

In the case of noisy data, incorporating e- congruence into this voting scheme is 
extremely difficult and complicated. Hence such a rigid bound on the point location 
error is normally not imposed for practical implementation. Clearly, if there is a large 
common substructure between two point sets P and Q, then in the transformation 
space there will be a large cluster of votes around the optimal transformation T. 
This is in contrast to the exact case, where all the votes corresponding to triplets of 
points belonging to LCP(P,Q) got accumulated precisely with the transformation 
T. Hence although the voting process is exactly the same as in the exact case, 
instead of trying to identify the transformation with the largest number of votes, 
the transformation space is searched for a large cluster of votes. Such a cluster is 
associated with the common substructure between P and Q. 


2.3 Geometric Hashing 

Given two point sets P and Q, the geometric hashing scheme divides the common 
point set identification problem into two phases, a preprocessing phase and a query 
phase. This technique is useful in cases where a point set P is compared with a 
large number of point sets Q u Q„, to identify LCP{P } Q t ) i s 1, 2, . . . , n. In 

such a situation, the preprocessing stage is done on P, independently of any of the 
QiS. The basic idea of this preprocessing is to store in a database, a representation 
of the point set P that is invariant under rigid transformation. By doing so, the 
representation of any Q { when computed during the query phase, will present some 
similarity with the transformation invariant representation of P. 


2.4 Related Work 
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There can be many ways in which this transformation invariant representation 
is calculated. We mention here a simple approach used in [PA94]. Since three non 
collinear points in space specify a plane and hence a reference frame, for every three 
points (pi,P 2 ,P 3 ) of P let R(p llP2 ,p 3 ) be a reference frame defined by these three 
points. For each such reference frame R(pi, P2 , P3 ), formed by triplets of points of P, 
the coordinates of all other points p € P are calculated in this reference frame and 
these coordinates are set as an index to a hash table where the couple (■ft(pi,p 2 ,p 3 ),p) 
is stored. 

During the query phase, as in the preprocessing phase, for every reference frame 
R(q ug2 ,q 3 ) (ft € Q, i = 1,2,3) the coordinates of all points q 6 Q, in this reference 
frame, are calculated. These coordinates are used to retrieve the pairs (R(p u p 2 , P3 ),p) 
from the hash table and for each such retrieved pair, a vote is casted for the triplet 
(Pi i ViiPz)- Clearly, all triplets having a large number of votes belong to LCP(P, Q). 
Randomization can again be used in the same way as was done in the case of 
alignment and voting. Firstly during the preprocessing phase, instead of using the 
entire point set P to calculate the reference frame, a randomly chosen subset RofP 
can be used, and the hash table be constructed. If any three points of LCP(P, Q ) 
belong to R then these points can be identified by the algorithm as they will receive 
a large number of votes. Moreover during the query phase also, all possible reference 
frames need not be considered. Depending on the sizes of the randomly sampled 
subsets, the running time decreases compared to the deterministic case. 


2.4 Related Work 

Apart from the approaches mentioned in the previous sections, algorithms based 
on clique detection have been studied widely in computational chemistry for this 
common substructure identification problem. For example DISCO [MBD + 93] builds 
a correspondence graph G from two molecules. The nodes of G are all possible atom 
pairs, and an edge is created if pairs of atoms can be matched simultaneously. A 
clique detection algorithm is then used to find cliques in G. These correspond to 
common substructures in the two molecules. Although maximum clique detection 
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is NP-hard [GJ80], this algorithm works well in practice [MBD+93, Wil95]. Similar 
graph theoretic techniques were proposed in [MARW89] for protein alignment, where 
it was suggested to build an undirected, labelled and fully connected graph, whose 
nodes correspond to the linear representation of the secondary structure elements 
of the proteins (i.e. helices and strands), and the edges correspond to the angles 
between these elements. For finding the common substructure between two proteins, 
a subgraph isomorphism algorithm is used on the corresponding graphs. 

The geometric hashing technique was used by Fisher et al. [FBNW92, FNW92] 
and also by Pennec at al. [PA94] for matching protein structures. Randomized 
versions of alignment were used by Finn et al. [FKL + 97] for identifying common 
substructures in drug molecules and by Akutsu [Aku96] for the protein structure 
alignment problem. Lesk, Pascarella and Argos, and Rao and Rossmann proposed 
iterative improvement methods [Les91, PA92, RR73]. Vriend and Sander developed 
a greedy method in which small fragments of protein molecules were assembled into 
larger structures [VS91]. Taylor and Orengo developed a double dynamic program- 
ming technique [T089], and §ali and Overington developed a stochastic method 
using probability density functions [§094]. However all of these methods have one 
or more limitations [HOS + 92] and most are not systematic, but based on some 
heuristic. Moreover as already mentioned, even though for practical purpose these 
methods may suffice in some cases, none, except the algorithms presented in refer- 
ence [Aku96] give a theoretical guarantee of the quality of the output. 


Chapter 3 


Approximation Algorithms for 3-D 
Common Point Set Identification 


In this chapter we first formalize the abstract geometric problem that was introduced 
in Chapter 1, and then give various approximation algorithms for the problem that 
arise out of this formalization. All of our algorithms are based on the alignment 
scheme described in the last chapter and have the underlying property that some 
aspect of the output has a guaranteed quality. There are two of these aspects that 
we consider here. The first is approximating the size of the largest common point 
set between two given point sets A and B, without violating the point location error 
e, and the second is finding out a common subset where e is satisfied only approxi- 
mately, both having guaranteed approximation ratios. As mentioned in Chapter 1, 
the second aspect was also addressed by Akutsu [Aku96], where, if S C A is the 
largest common point set between A and B under e-congruence, then his algorithm 
outputs a subset S' C A, with J-S'l > \S\ and S' being 8e-congruent to some subset of 
B. Reducing this factor of 8 involves a large constant factor in the time complexity 
of the algorithm. We reduce this factor to 2 in our algorithm without incurring any 
such large constant factor. The approximation approach that we take here extends 
the ideas from [GM094, Aku96] and also from [Sch92, Sch88] where approximate 
decision algorithms for approximate congruence of two equal cardinality point sets 
were presented. 
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3.1 Statement of the Problem 

Definition 3.1 A point set S is e-congruent to a point set S' if there exists an 
isometry I and a bijective mapping l : S -> S' such that for each point s € S, 
d(I(s),l(s)) < e, where d(- } •) is the underlying metric. 


Definition 3.2 Given point sets A , B and real numbers e > 0 and 0 < a < 
1, a-LCP(A, B, e) is a subset S of A with \S\ > amin(|A|, \B\) such that 5 is 
e-congruent to a subset of B. 


This version of LCP was introduced by Finn et al. [FKL + 97j, where is was refered 
to as LCP-a, with the motivation that this captures the primary application quite 
well, since in practice the common substructure between two molecules must only 
be large enough. Clearly, for any e, there exists a a max (e) such that a-LCP(A, B, e) 
exists for all a < a max (e) and for any a > cx Tnax (e), a- LCP {A, B, e) does not exist. 
Similarly, for any 0 < a < 1, there exists e m i n (a) such that a-LCP(A, B, e) exists 
for all e > e mi „(a) and for any e < e min (a), a-LCP(A, B, e) does not exist. Consid- 
ering this, we address the following problem : 

Input : 3-D point sets A, B and a real number e > 0 
Output : a ma x(e)-LCP(A,B 1 e) 

Unless otherwise mentioned, from how onwards a point set refers to a point 
set in 3-D, and any isometric transformation is a composition of just rotation and 
translation, not including mirror image. This restricted definition of an isometry 
does not result in any loss of generality, because isometry including mirror image 
just increases the computation time of any of our algorithm by only a constant 
factor. For any two point sets A and B, now the algorithm has to be run once with 
A and B and then with A and a mirror image of B on a plane that r an be chosen 
arbitrarily [Aku96]. 
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3.2 Approximating the Size of the Largest Com- 
mon Point Set 

In this section we state an algorithm which, when given point sets A, B and a real 
number e > 0, outputs a point set ai-LCP(A, B, e) and a real number a u such that 
011 5s &max(,€) ^ 

For point sets A and B, let X be the isometry and l be the bijective map- 
ping corresponding to Oi-LCP(A,B,e min {aj) Q A. Let pi be any point of a- 
LCP(A,B,e m in(a)). Let p 2 be the point belonging to a-LCP(A, JB, e m tn(a:)) which 
is farthest from pi and pz be the point of a-LCP(A,B,e m i n (a)), such that the 
perpendicular distance from pz to the line passing through p\ and pz is maximized. 

Now, let Ti be the translation that takes the point pi to R\ be the rotation 
about the point Ti(pi), that causes Ti(pi), Ti(p 2 ) and Z(p 2 ) to become collinear, and 
J R 2 be the rotation about the Ri (T x (pi))Ri (7\ (p 2 ) ) axis, that causes J?i(Ti(pi)), 
i?i(Ti(p 2 )), R\ (Ti(pz)) and l(pz) to become coplanar. Let X' be the isometric trans- 
formation which is the composition of T l5 Ri and i? 2 , i.e. X'(p) = ,R 2 (jRi(Ti(p))). 
Then the following lemma follows immediately from [GM094]. 

Lemma 3.1 The isometry X' and the bijective mapping l correspond to 
a-LCP(A , B , 8e m i n (a)) where T , l, A, B, a are as described above. 

Proof: If the set X(A) is translated such that the point pi becomes coinci- 

dent with l(pi), then each point p € X (ot-LCP(A,B,e min (a))) is within a 'distance 
2c mi „(o;) from its corresponding point l(p) € B. Now if the resulting set is ro- 
tated about the point pi such that p x , p 2 and l(p 2 ) become collinear, then the point 
P 2 moves by atmost a distance of 2e m i n (a). Since p 2 is the farthest point of oc- 
LCP(A , B, € min (a)) from the centre of rotation, any point of a-LCP{A, B, e mi „(a)) 
will move by a distance of atmost 2e min (a). Finally, when the resulting set is ro- 
tated about the pipj a>ds such that p\, P 2 , Pz l(pz) become coplanar, the point 
Pz moves by a distance of atmost 4e m tn(&). Since pz is the farthest point of oc- 
LCP(A, B, e m i n {oi)) from the pi pi axis, any other point of a-LCP(A, B, e m in(a)) 
moves by a distance of atmost 4e m i n (a). So, any point of a-LCP(A, B, e m i n (a)) was 
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moved by atmost £ m j n (a) by translation, 2e m i n (oc) by the first rotation and 4e m j n (a) 
by the second rotation. Since this point was originally atmost a distance of € m i n (o/) 
away from Z(p), now it is atmost a distance of 8e m i n (oi) away from l(p). This proves 
the lemma. I 

Definition 3.3 If P is a triplet of points and Q is another triplet of points, let T PQ 
be the isometric transformation, as described above, i.e. if P — (p»i»P* 2 ,P* s ) and 
Q = then T PQiPh) = Qhi TpqM, T PQ (jp i% ). and q h are collinear and 

T PQ (p h ), T PQ (p i2 ), T PQ (p i3 ) and q h are coplanar. 

D efini tion 3.4 For point sets A, B, isometry X : A — + B and a real e > 0, let 
G(X, e, A, B) be a bipartite graph (U U V, E ) where U and V represent the points of 
A and B respectively and if u € U is the node corresponding to a € A and v € V 
corresponds to 6 € B, then E = {(«,u) | d(l(a),b ) < e}. 

Next, making use of the isometry stated in Lemma 3.1, we first give an approxi- 
mate decision algorithm about the existance of a -LCP(A, B , e) where point sets A , 
B, and real numbers a and e are input to the algorithm. This decision algorithm 
is called approximate , because it can take a decision only for values of (a, e) for 
which e is not too close to e min (a;). When e is too close to e m i n (a), the algorithm 
may return DON’T KNOW and such values of € are said to constitute the inde- 
cision interval. However, whenever the algorithm returns YES or NO, the answer 
is correct. This algorithm has an indecision interval equal to [|c m i„(a),8e min (a)). 
Making use of this decision algorithm, we construct an algorithm which approxi- 
mates oc max (e)-LCP(A, B, e) in the way mentioned at the begining of this section. 
We then calculate the running time of this algorithm and finally show that the qual- 
ity or accuracy of this approximation is related to the size of the indecision interval 

[ g €»wn (^ » 8^min ) • 


Algorithm 3.X 

Input : Point sets A, B, real numbers c > 0, 0 < a < 1. 

Output : If a-LCP{A,B,e ) exists, then YES or DON’T KNOW, else NO or 
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DON’T KNOW. 

for all triplets P = (Pi n Pi 2 ,Pi 3 ) from A 
for all triplets Q = (g^, qj 3 ) from B 

{ 

compute Tpq\ 
construct G(Tpq , e, A, B); 

if G(!Tpq, e, A, B) has a matching of size greater than or equal to a min(|«A[, \B\) 
then return YES; 

} 

decision = NO; 

for all triplets P = (pn,Pi 2 ,Pi 3 ) from A 
for all triplets Q = (q^, q j2 , q j3 ) from B 

{ 

compute Tpq‘, 
construct G(Tpq, 8e, A, B ); 

if G(Tpq , 8e, A, B ) has a matching of size greater than or equal to a min(|A|, |B|) 
then decision = YES; 

} 

if (decision == NO) then return NO else return DON'T KNOW; 

Theorem 3.2 Algorithm 3.1 always returns the correct answer about the existance 
of a-LCP(A,B,e) if 

(i) e > 8 € mi „(a), or 

(ii) e < fe m i„(aO 

and it either returns the correct answer or returns DON’T KNOW if g£mm(oO ^ 
e < 8e min (a). 

Proof: If an isometry T and bijective mapping l correspond to a-LCP(A, B, e), 

then I and l correspond to any a-LCP{A, B, <f) where e' > e. From Lemma 3.1, it 
follows that Algorithm 3.1 finds an isometry X and a mappinf Z which enables in 
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finding an a-LCP(A, B, 8e min (a)). Since X and l also correspond to a-LCP(A, B, e) 
for any e > 8e min (a), the algorithm returns YES for all € > 8e m i n (a). 

For any e, the algorithm returns YES iff G(Tpq, e, A, B ) has a matching of size 
greater than or equal to amin(|A|, |B|). Hence all YES answers are correct. 

For any e < |e min (a), no isometry enables in finding a-LCP(A,B,8c). Hence 
for any e < |e min (a;), the algorithm always returns NO. 

Finally, since the algorithm finds an isometry and bijective mapping correspond- 
ing to a-LCP(A, B, 8e min (a)), it can return NO only if e < e min (a). Hence all NO 
answers are also correct. | 

Lemma 3.3 If \A\ = m and |S| = n, then Algorithm 3.1 has time complexity 
0(m 4 n 4 v / m + n). 

Proof: There are 0(m 3 ) triplets of A and 0(n 3 ) triplets of B. Constructing the 

graph G(T P Q,e,A,B) takes time 0(mn). Computing the maximum matching in 

G(T P q, e, A, B) can be done using Hopcroft and Karp’s algorithm [HK73] in time 
0(mny/m + n). | 

Algorithm 3.2 

Input : Point sets A, B and a real number e > 0. 

Output : Real numbers 0 < og < a u < 1. 

M = 0; 

for all triplets P = (p h ,p i2 ,p h ) from A 
for all triplets Q = (q h ,q h ,q j3 ) from B 
{ 

compute Tpq-, 
construct G(T PQ , e, A, B); 

M' = size of the maximum matching in G(T PQ ,e, A, B); 
if [M 1 > M) then {M = M'\ T — T PQ ; } 

} 

ai = Af/min(lA|, |J3|); 
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M = 0; 

for all triplets P = ( Ph,Pi 2 ,Pi 3 ) from A 
for all triplets Q = (qj^qj^q^) from B 

{ 

compute Tpq ; 
construct G(Tpq, 8e, A, B ); 

M' = size of the maximum matching in G(TpQ,8e,A,B ); 
if (M 1 > M ) then M = M'\ 

} 

ot u = (M + l)/min(|A|,|B|); 
return (o>i } a u ); 

This algorithm is based on the same idea as that of the decision Algorithm 3.1. 
Clearly, for a fixed e, a L is the largest value of a for which Algorithm 3.1 would 
return YES and a u is the smallest value of a for which Algorithm 3.1 would return 
NO. 

Theorem 3.4 If \A\ - m and \B\ = n, then Algorithm 3.2 has time complexity 
O (m 4 n 4 s/m + n) . 

• Proof: Same as that for Lemma 3.3. | 

Theorem 3.5 Given point sets A, B and a real number e > 0, Algorithm 3.2 returns 
real numbers 0 < on < a u < l, such that on < a max (e ) < a u and 

(i) on > max {a : e > 8e m j n (o:)} 

(ii) a u < min {a : e < |e min (a)} 

Proof: Clearly, Algorithm 3.1 returns YES for a = aj and NO for a = a u . 

Hence the first part of the theorem is obvious, since on-LCP(A,B,e) exists and 
a u -LCP(A,B,e) does not exist, along with the fact that, for all a < a max (c), a- 
LCP{A, B, c) exists and for any a > oc max {e), a-LCP(A, B, e) does not exist. 

The second part of the theorem follows from Theorem 3.2 along with the fact 
that if a x < a 2 then, e m m{o!i) < t m in{oL 2 ). I 
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3.3 An Exact Algorithm for Finding the Largest 
Common Point Set under Rotation 

Now we describe an algorithm which on given point sets A, £, a real number e > 0, 
and a fixed point p, finds a max {e)-LCP{A, B, e) where the underlying isometry con- 
sists of only pure rotation about the point p. We shall make use of this algorithm 
in the next section to decrease the indecision interval of Algorithm 3.1 and thereby 
obtain a better approximation of OL max {e)-LCP{A, B, e.) under general isometry, com- 
pared to what is obtained by Algorithm 3.2. 

To understand this algorithm, consider e-balls around each point of the set B. 
As the set A is rotated about the point p, points of A move into and out of the e-balls 
of B. The problem is essentially that of finding the rotation for which maximum 
number of points of A are within distinct balls of B. If A and B were 2-D point 
sets, then a straightforward approach could have been as follows : for each point 
€ A and each point bj 6 B, let % be the angular interval for which a* is within 
the e-ball around bj, when the set A is rotated about the point p. The end points 
of these intervals are sorted to get a new set of intervals. Now, we move through 
this new set and for each angular interval in this set, construct the bipartite graph 
as was done in Algorithms 3.1 and 3.2, in which each edge (a,, bj) indicates that a* 
is within the e-ball around bj for this angular interval. The interval for which this 
graph has the largest maximum matching, results in finding the largest common 
point set between A and B under rotation about p. For 3-D point sets, we do not 
obtain angular intervals, but rather solid angles. So moving through these angles is 
not as straightforward as in the case of two dimensions. 

To describe the algorithm we make use of the following notation. If p is the fixed 
point about which the point set A is rotated, then for any point a, € A, S ai denotes 
the sphere of radius d(p, a,) centered at the point p, and for any point bj € B, S bj 
denotes the sphere of radius e centered at bj, i.e. denotes the e-ball around the 
point bj. S p will be used to denote a sphere centered at p. The radius of this sphere 
is less than the distance of p from the nearest point of either A or B. This implies 
that no point of either A or £ is within the sphere S p . Although this condition 
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is not necessary for the algorithm to work, it will make the situation simpler by 
reducing the number of different possible cases. If spheres S a . and S bj intersect, 
the intersection of their surfaces will be a circle [Lei43], which we denote by Cy. 
Consider the solid cone having p as its vertex and the circle Cy- as its base. This 
solid cone intersects with the sphere S p , to form a circular figure called a dome , on 
the surface of S p . This dome is indicative of the solid angle corresponding to which 
the point a { is within the e-ball around the point bj. Now consider a fixed point 
p' on the surface of S p . If the same set of rotations which are to be applied to the 
point a,{ for it to lie within the e-ball around bj, are applied to the point p 1 , then it 
traces out an exactly similar dome on S p . Let us denote this dome traced out by 
the point p' , by Dy (see Figure 1). 



Figure 1: Dome Dy resulting from points a* and bj 

If there is a rotation R of the set A, about the point p, such that d(R(ai),bj) < e, 
then obviously Dy # 0. Clearly, for any point on the dome Dy, there exists a 
rotation R such that d(R(ai),bj ) < e. From now onwards all our references to 
domes are to the circular figures on the surface of S p , traced out by the point p' as 
the point set A is rotated so that a point of A is within the e-ball of another point of 
B. Now consider every point on the surface of the sphere S p to be associated with a 
membership vector , which is indicative of all the domes to which the point belongs. 
Then each dome partitions the surface of S p into two regions, and all the domes 
arising out of the points of A and B define a partition of the surface of S p into a 
number of distinct regions, where a region is defined as a set of points having the same 
membership vector. Therefore, for any point on the region Dyy nDyy n. . .flDy k , 
there is a rotation R such that d(R(ai i ),bj l ) < e, l = 1, 2, This gives rise 
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to a bipartite graph, in which the nodes correspond to the points of A and £, and 
the edges consist of all pairs (a in b jt ), l = 1 , 2 ,..., A:. The region in which this 
graph has the largest maximum matching is our required region, because a rotation 
corresponding to any point in this region finds the largest common point set between 
A and B. To find the largest maximum matching, it is required to traverse through 
all the regions and find the maximum matching in the bipartite graph arising in 
each region. To systematically traverse through all these regions we make use of a 
sweep approach. Towards this, we sweep a plane h(t) : x = t, through the sphere S p , 
starting from its leftmost end and ending in its rightmost end. At any position x = t, 
this plane intersects with the sphere S p resulting in a circle C(t). This situation is 
illustrated in Figure 2. 



Figure 2: Sweeping the plane h(t) :x = t through the sphere S p 

Any such circle C(t), being a section of the surface of S p , intersects with some 
of the domes on the surface of S p . For any such dome Py, let 0y be the angle 
subtended by the intersection points of Py with C{t), at the centre of the circle 
C(t) (Figure 3). 

For all such intersecting Py, we compute the corresponding angular intervals 
%• We next sort the end points of these angular intervals, which as a result 
gives rise to a new set of intervals. If an interval is contained within the domes 


3-3 An Exact Algorithm for Finding the Largest Common Point Set under Rotation 25 



Figure 3: Angle % resulting out of the intersection of the dome Da with the circle 

C(t) 

Di^, Di Vl , . . . , An.;*) then it results in the bipartite graph corresponding to the 
region formed by the intersection of these domes. Hence if we find the maximum 
matching in all the graphs corresponding to each interval for this position of the 
plane, and repeat this process for all positions of the plane by sweeping it through 
S p , then the largest maximum matching over all possible graphs is our required 
result. However, as is normal in any sweep algorithm, this (impossible) task of 
generating an infinite set of planes can be avoided. Let x x ,X 2 ,...,x n be the sorted 
^-coordinates of all the left and right end points, and all possible intersection points, 
of all the domes Dy. Then for any position of the plane between two consecutive XjS, 
i.e. h(t) : x — t, t£ (x,, Xi +X ) (i = 1, 2, . . . , n— 1), the same set of graphs arise out of 
the intersection of the circle C(t) with the domes lying on it. Hence it is sufficient to 
consider only one position of the plane between any two of these successive XjS. If A 
and B are point sets of cardinality m and n, then there are 0{mn ) domes on S p and 
hence 0(m 2 n 2 ) intersection points. For each position of the plane h(t) : x = t, there 
are 0(mn ) angular intervals in the corresponding circle C(t). Hence this scheme 
needs 0(m 3 n 3 ) invocations of the graph matching algorithm. However, note that 
0(mn) domes divide the surface of the sphere S p into 0{m 2 n 2 ) regions. Hence the 
graph matching algorithm should be invoked only 0{m 2 n 2 ) times. For this we note 
that if S x . S 2 , ■ • • ,Sk axe sets of domes, such that for each Si, there is a region n 
which is contained in all the domes of Si, then it is sufficient to do a graph matching 
only for the regions r,-, i = 1, 2, . . . , k. Obviously the sets Si, S 2 , • . . , S* need not 
be disjoint. Figure 4 illustrates this point. It is so, because, say if D tlJ1 and Dj 2J - 2 
intersect, then the graph arising out of the common region will have its edge set as 
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{(0^,6*), (a* 2 , 6 j 2 ) } , but the other two graphs will have their edge sets as {(<^,6*)} 
and {(a i2 , b h )} only. Hence the size of the maximum matching in the graph arising 
out of the common region will be atleast as large as the maximum matching in 

S,- (C fCj, C 4 ) 

Sj « (C,.C,C S > 

S } . IC4.C s l 


Figure 4: It is sufficient to do graph matching corresponding to rj, r 2 , r 3 and r 4 only 

Therefore if A U1 , A 2 j 2 » • • • , A k j k be a group of domes such that there is a region 
r which is contained in all of them, our objective is to identify this region and all 
the domes containing it, using our sweep algorithm. Observe that there are three 
possible cases : i) The region r does not share its boundary with atleast one of 
the domes i.e. r is in the interior of the dome (Figure 5 (a)). ii) One of the domes 
constitute the region r (Figure 5 (b)). iii) r shares its boundary with all the domes 
(Figure 5 (c)). 

© 

(a) (b) (c) 

Figure 5 : Three possible arrangements of domes 

Let us consider the third case first. If pi,p 2 ,.-.,Pn be the intersection points 
lying on the boundary of r, then clearly all of these points have their membership 
vectors equal equal to that of any point in the region r. Hence it is sufficient to do a 




graphs corresponding to the two other regions. 
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graph matching corresponding to any of these points. The same scheme works in the 
first case also, if we consider all possible intersection points. This is because all the 
intersection points lying on the boundary of r also lies inside the dome which does 
not share its boundary with r. However, this scheme of doing the graph matching 
at the intersection points fails in the second case. To take care of this, note that any 
point lying on the boundary of the dome which constitutes the region, can fulfill our 
requirements. For any such point, we can find out all the domes to which this point 
belongs and construct the corresponding graph, which will be the graph underlying 
r. We then do a graph matching corresponding to this graph. 

Now we give the overall algorithm. Observe that the membership vector of the 
sweep plane indicating the domes which are intersected by the plane, changes only in 
two situations : (i) the sweep-plane just crossing the left end-point of the dome, after 
which this dome starts intersecting with the sweep-plane (ii) the sweep-plane just 
crossing the right end-point of the dome, after which this dome ceases to intersect 
with the sweep-plane. Our event point schedule , i.e. is the sequence of abscissae 
ordered from left to right, which define the halting positions of the sweep-plane, are 
made up of the ^-coordinates of the left and right end points, and all the intersection 
points of all the domes lying on S p . 

When the event point is the left end-point of a dome, we update the membership 
vector of the sweep-plane to indicate that this dome now intersects with the sweep- 
plane. Next we construct the bipartite graph corresponding to this point making 
use of the information in the sweep-plane membership vector, because this point can 
lie only on the subset of all the domes which intersect with the sweep-plane. If the 
event point is contained in the domes Aui, Asia* • > A k j k , then the nodes of the 
bipartite graph correspond to the points of A and B and the edges join pairs a € A, 
bj t € B, l = l,2....,fc. We then compute the size of the maximum matching in 
this graph. If the event point is an intersection point of a number of domes, then 
we construct the bipartite graph corresponding to this point as was done in the 
previous case, and find the maximum matching in this graph. Finally, if the event 
point is the right end-point of a dome, then we just update the membership vector 
of the sweep-plane to indicate that from now this dome ceases to intersect with the 
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sweep-plane. The largest maximum matching obtained over all the graphs is our 
required result. Now we formally state our algorithm for finding the largest common 
point set under pure rotation, and the consequent theorem concerning its running 
time. 

Algorithm 3.3 

Input : Point sets A, B, a fixed point p, and a real number e > 0. 

Output : ctmax (c) for pure rotation of the set A about the point p. 

for each Oj € A 
for each b{ E B 
compute Dij ; 

^-structure = sorted ^-coordinates of the left and right end-points and intersection 
points of all the domes Dy ; 

M = 0; 

for each t € rr-structure 

{ 

if (t = = left end-point of a dome) then 

{ 

update sweep-plane membership vector; 

construct graph G corresponding to the point t using information in the sweep-plane 
membership vector; 

M' = size of the maximum matching in G; 

} 

if (t == intersection point of domes) then 

{ 

construct graph G corresponding to the point t using information in the sweep-plane 
membership vector; 

M r = size of the maximum matching in (7; 

} 

if (t == right end-point of a dome) then 
update sweep-plane membership vector; 
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if (M' > M) then M = M'\ 

} 

return M/min(|A|, |J3|); 


Theorem 3.6 If |A| = m and | J5 1 = n, then Algorithm 3.3 has time complexity 
0{m z n z y/m+ll) ■ 

Proof: The m points of the set A and n points of the set B result in 0(mn ) 

domes on the surface of S p . Hence there are 0(mn) end-points and 0(m 2 n 2 ) in- 
tersection points. Corresponding to each of the 0(m 2 n 2 ) event points, constructing 
the bipartite graph takes 0(mn) time and finding the maximum matching takes 
0(mnVm + n) time. Hence the overall time complexity is 0{m z n z y/m + n). | 


3.4 Reducing the Indecision Interval 

It was shown in Section 3.2 that the indecision interval of the decision Algorithm 3.1 
is [|€ m i n (a), 8e m i n (a)). In this section we shall make use of Algorithm 3.3 and reduce 
the indecision interval to [\e m i n (a), 2 e m i n (<*))• Using this new decision algorithm, we 
next give an algorithm which better approximates the size of cx max (e)-LCP(A, B, e) 
compared to Algorithm 3.2. For this we first give a lemma which is based on the 
same idea as that of Lemma 5 of [Sch92], but adapted for our particular problem. 
Here, for arbitrary points a and b in space, we use t a b to denote the translation that 
maps a to b. 

Lemma 3.7 Let isometry X, which is a composition of translation and rotation, 
and a bijective mapping l, correspond to a-LCP(A, B,e). Let a € a-LCP(A,B,e). 
There exists a rotation R about the point X(a), such that R and l correspond to 
a~LC P(t a x(a){A), B , e). Let b be an arbitrary point in space. There is a rotation R! 
about the point b such that R! and l correspond to a-LCP(t a b(A), B,e + d(b,X(a)). 

Proof: The first part of the lemma is obvious. From this it follows that X = 

R°tax(a) (where o is used to denote composition carried out from right to left). For 
the second part of the lemma, we show that R' = tj(a)6 ° R ° t~ l x{a)b- 
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For any point p €• a-LCP(A, B , e), 


d(R'ot ab (j>),l(p)) = 


< 

< 


d( tj( a )b ° R°t 1 1(o)6 ° t ab (p), tfp) ) 
d(tz(a)b 0 R° tal(a) {p ) > Kp) ) 
d{tx(a)b°Z(p), Kp) ) 
d(l(p),Z{p)) + d(l(p),h( a )b°Z(p)) 
c + d(b,Z(a)) 


I 

Algorithm 3.4 

Input : Point sets A, B, real numbers c > 0, 0 < a < 1. 

Output : // a-LCP(A, B, e) exists, then YES or DON'T KNOW , else NO or 
DON’T KNOW. 

for each point a G A 
for each point b € B 

{ 

a ' = Algorithm 3.3(t ab (A),B,b,e)] 
if (o' > O') then return YES; 

} 

decision = NO; 
for each point a € A 
for each point 6 E B 

{ 

a' = Algorithm 3.3(t ab (A),B,b,2e); 
if (a' > a) then decision = YES; 

} 

if (decision = = NO) then return NO else return DON’T KNOW; 


Theorem 3.8 / A and B are point sets of cardinality m and n, then Algorithm 8*4 
runs in time 0 (m 4 n 4 \/m + n) and always returns the correct answer about the ex- 
istence of a-LCP(A,B,e) if 
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(i) e > 2e min (a), or 

(H) 6 < 2^tnin(.®) 

and it either returns the correct answer or returns DON’T KNOW if le min (a) < 
C < 2e m i n (a). 

Proof: There are 0(mn) invocations to Algorithm 3.3. Hence the time complex- 

ity of Algorithm 3.4 is 0(m 4 n 4 \/yn"+"n). Let isometry X and a bijective mapping 
l, correspond to a-LCP(A,B, e min (a)). Let a € a-LCP(A, B, e min (a)). It follows 
from Lemma 3.7 that there exists a rotation R about the point 1(a) , such that R and 
l correspond to a-LCF’(t al ( a )(A),B,e min (a) +d(l(a),X(a))). Since d(l(a),I(a)) < 
tmin (ot ) , R and l correspond to a-LCP(t a ^ a )(A),B, 2e min (a)). The rest of the proof 
is similar to that for Theorem 3.2. | 

Now in the same way as was done in Section 3.2, this decision algorithm can be 
modified to give an approximation algorithm, which when given point sets A, B, 
and real number e > 0, outputs on and a u such that at < a max < a u , and finds the 
common point set at-LCP(A,B,e). In contrast to Theorem 3.5, in this case at > 
max {a: : 2e min (a)}, and a u < min {ct : e < ^Cmin (<*)}• Hence this gives a signif- 
icantly better approximation of a max (e)-LCP(A,B,e). Clearly this approximation 
algorithm will also have a time complexity equal to that of the decision algorithm. 


3.5 Approximately Satisfying the e-Constraint 

In all the algorithms presented so fax, the constraint imposed by the point location 
error e was strictly satisfied. In this section we present an algorithm, which on given 
point sets A, B, and real number e > 0, instead of approximating a max (e) as was 
done till now, approximately satisfies the constraint imposed by e. Towards this, the 
algorithm outputs a subset S C A which is 2e-congruent to some subset of B, and 
\S[ > \a max (e)-LCP(A, B, e)|. This approach was also followed by Akutsu [Aku96] 
and his algorithm outputs a subset which is 8e-congruent to some subset of B. We 
have made use of our exact algorithm for finding a max (e)-LCP(A, B, e) under pure 
rotation, to improve this bound to 2e. 
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Algorithm 3.5 

Input : Point sets A, B, real number e > 0. 

Output : A real number oc. 

a = 0; 

for each point a & A 
for each point b e B 

{ 

a! = Algorithm 3.3(t a b(A), B,b,2e); 
if (o' > a) then a = a'; 

} 

return a; 

Theorem 3.9 Given point sets A, B, and a real number e > 0, if Algorithm S.S 
returns a real number a, then there exists a subset S of A with jS| *= amin(\A\, |£|), 
which is 2e-congruent to some subset of B and a > £*»»«*(€). 

Proof: Clearly, it follows from the construction of Algorithms 3.3 and 3.5 that if 

Algorithm 3.5 returns a, then there exists a subset S of A with |5| » amm(JA|, |£|). 
Let isometry X and bijection mapping l correspond to a m<tx {a)-LCP{A, B, «). If a € 
a max (e)-LCP(A,B,e), then from Lemma 3.7 it follows that there exists a rotation 
R about the point X(a) such that I - R o tax( a ), and the isometry T a# R o $**{*) 
correspond to a-LCP(A,B,2e), where a > amax(e). Since Algorithm 3.5 loops 
through all possible pairs of points of A and B, certainly Algorithm 3.3 gets called 
for a pair a, 1(a) where a 6 a max (€)-LCL(A, B,e). Hence the theorem holds. I 


Theore m 3.10 If |A| = m and |£| = n, then Algorithm S.S has time complexity 
0 (m 4 n 4 ^J(m + n). 


Proof: There are mn calls to Algorithm 3.3. | 
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3.6 Running Time Versus Size of the Indecision 
Interval 

The algorithm that we gave in Section 3.2 had an indecision interval of [fe min (a:), 8e min (a 
Using our exact algorithm of Section 3.3 we reduced this indecision interval to 
[|e m tn(a;), 2e m i n (o;)) in Section 3.4, and showed that it leads to a better approxi- 
mation of Oi max (e). Making use of a technique first described in [Sch92], in this 
section we show that this indecision interval can be made arbitrarily small, however, 
at the expense of increased running time. Given point sets A, B, and a real number 
e, we can reduce the indecision interval to [e min (a:) — j,e min (a) + 7), by slightly 
modifying Algorithm 3.4. But doing this introduces a term (e/7) 3 in the running 
time of the algorithm. In this algorithm we cover the e-balls around each point of 
the set B with balls of radius 7. Let B' be the set of points which are the centers 
of these 7-balls. 

Algorithm 3.6 

Input : Point sets A, B, real numbers e>0, 0 <a<lj 0<7 <e. 

Output : If a-LCP(A,B,e) exists, then YES or DON’T KNOW, else NO or 
DON’T KNOW. 

for each point a € A 
for each point b' € B' 

{ 

a! = Algorithm 3.3(toy (A),B, b', e); 
if (of > oc) then return YES; 

} 

decision = NO; 
for each point a £ A 
for each point V £ B' 

{ 

o' = Algorithm 3.3(<«*(A), B, b’, e + 7); 
if (a' > a) then decision = YES; 
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} 

if (decision = = NO) then return NO else return DON'T KNOW; 


Theorem 3.11 Algorithm 3.6 always returns the correct answer about the existence 
of a-LCP(A, B, e) if 

(i) £ > £min (o) 4" 1 ) 

(ii) € < £min (®) 7 

and it either returns the correct answer or returns DON’T KNOW if e min (a) — 7 < 
e < € m i„(a) +7. 

Proof: Let isometry 1 and bijective mapping l, correspond to a~LCP(A t B, t m i n (a)). 

Let a € a-LCP(A, B, e m ,„(a)). Then clearly the isometry X maps o into the c-ball 
around 1(a), and hence a 7 -ball inside this e-ball. Let the center of this 7 -ball be 
V. Prom Lemma 3.7 it follows that there exists a rotation R, which along with l 
correspond to a-LCP(tat/(A),B , e m j„(o!) + 7 ). The rest of the proof is exactly along 
similar lines as that for Theorem 3 . 2 . | 

Theorem 3.12 If \A\ = m, |J3| = n, then Algorithm 3.6 has a time complexity 
0((e/ 7 ) 3 m 4 n 4 v'm + n) . 

Proof: ( 2 e/ 7) 3 balls of radius 7 are sufficient to cover each e-bali around the 

points of B. So for each point of B, there are 0((e/ 7 )®) additional iterations com- 
pared to that in Algorithm 3 . 4 . | 

In the same way, as was mentioned in Section 3.4, this unproved decision algo- 
rithm leads to an even better approximation of a max (e). However the running time 
of such an algorithm increases by 0 ((e/ 7 ) 3 ). By making 7 arbitrarily small, we can 
now get ai and a u which will be closer to <Xmax(e) compared to what was obtained 
in Section 3.4. 

This same technique can be applied to Algorithm 3.5 to reduce the factor of 2 
to any 8 > 1 , by appropriately choosing 7 . But here also we shall get a factor of 
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(e/ 7) 3 in the running time. This is however significantly better than the algorithm 
proposed by Akutsu [Aku96] which obtains the same result but introduces a factor 
of (e/ 7) 9 in the running time. 



Chapter 4 


Improvements Using Geometry, 
Randomization, and Structure of 
the Molecules 


In this chapter we present three different modifications of the basic algorithms of 
the last chapter, which improve their running time. The first two are completely 
general, while the third relies on some structural properties of protein molecules. 


4.1 Using an Approximation Algorithm for Max- 
imum Matching in Bipartite Graphs 

All the algorithms presented in the last chapter make use of the Hopcroft and 
Karp’s algorithm [HK73] for finding the maximum matching in a bipartite graph. 
Clearly, any improved algorithm for maximum matching would result in improving 
the overall time complexity of our algorithms. It was shown by Efrat and Itai [EI96] 
that if the nodes of a bipartite graph are points in some d-dimensional space, and 
the edges are pairs of points which are within some specified distance of each other, 
as in our case, then atleast under certain circumstances, algorithms with better 
time complexities can be obtained. This was done by making use of the geometric 



4.1 Using an Approximate Algorithm for Maximum Matching in Bipartite Graphs 37 


nature of the problem, rather than trying to solve it from a purely graph-theoretic 
perspective. Work related to this was also done by Vaidya [Vai89j, and by Agarwai, 
Efrat and Sharir [AES95). 

In this section we show that making use of the scheme of Efrat and Has [EI96J 
along with a data structure proposed by Arya et al. [AMN+94] for answering nearest 
neighbor queries for points in 5R d , we can reduce the time complexity of all our 
algorithms. However, the resulting algorithms have an increased indecision interval 
and the point location error e increases by a constant factor. 

4.1.1 Maximum Matching Using Geometry 

For the sake of completeness, in this subsection we outline parts of Efrat and Itai’s 
scheme that we shall make use of in our algorithm in the next subsection. Given 
point sets A, B, and a real number e > 0, consider the graph G = (A U B, E) where 
E consists of all pairs (a, b) a € A, b € B, for which d{a, b ) < €. Our problem is to 
find a maximum graph-matching in G. 

First we recapitulate some basic definitions. A matching M of G = (A U B, E) 
is a subset of E such that no vertex of G is incident on more than one edge of M. 
The vertices incident on M are called matched and the remaining are referred to as 
exposed. A path tt = (vi, v 2 , . . . , v 2 t) is called an alternating path if t>i is an exposed 
vertex of A, (va,ifc+i) e M and 6 E \ M (* = 1,2, . . . >*). Obviously 

the odd vertices of 7r belong to A and the even vertices to S. This path is called 
an augmenting path if v 2 t is an exposed vertex. If ir is an augmenting path, then 
M' = (M\ tt) U (t r\M) is also a matching and \M’\ a* \M\ + 1. Moreover, a theorem 
due to Berge [Ber57] states that a matching is maximum if and only if there is no 
augmenting path. Hence an algorithm for maximum matching can start with an 
empty matching and augment it by finding augmenting paths, until none is found. 

Given a matching AT, Efrat and Itai’s algorithm finds all shortest augmenting 
paths by first conducting a breadth-first search over the vertices of G to construct a 
layer graph C consisting of layers L\ t L 2 , . . . , L 2 t- The first layer Lu contains all the 
exposed vertices of A; L 2i contains all vertices of B not appearing in Uj< 2 * and 
connected through edges of G to some vertex of L 2i -i. If L 2i contains any exposed 
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vertices, then this is the last layer. Otherwise L 2i+ 1 is constructed, which contains 
all vertices connected through the matching M to vertices in L 2i . The layer graph L 
consists of the vertex set edges of M that connect vertices of L 2j to vertices 

of L 2j+ 1 , and edges of G that connect vertices of L 2j -i to vertices of L 2j . 

Towards finding a maximum matching, an algorithm due to Dinitz [Din70] finds 
a maximal set of edge-disjoint augmenting paths by conducting a depth-first search 
of the layer graph C. Efrat and Itai’s algorithm takes advantage of the geometric 
setting of the graph to improve the efficiency of Dinitz’s algorithm. For this an 
abstract data structure D, (5) is used, for a set of points S. Using this data structure 
the set of vertices of each layer L { is found, without explicitly constructing all the 
edges of £. Given a set of points S and a real number e > 0, the data structure 
should support two operations : 

(i) neighbor £ (S, q) : For query point q , return a point s € S whose distance from q 
is atmost e. If no such s exists, then neighbor (S, q) = 0. 

(ii) deleters, s) : Delete the point s from S. 

Let T(|5j) denote an upper bound on the time of performing one of these operations. 
At this moment let ns disregard the time required to construct this data structure, 
since it will be shown in the next subsection that this time in our problem is bounded 
by 0(nT(n)), where S is a point set in 3-D with cardinality n. Hence it will not 
influence the overall time complexity. Using this data structure, Algorithm 4.1 shows 
how the layer graph C is constructed. It is to be noted that each matched vertex 
of A is reached in constant time from its pair in M, and each vertex of B is found 
atmost once by a query of neighbor^!), •) and also deleted from D atmost once. 
Hence constructing C requires 0(nT(n )) time. 

Algorithm 4.1 (Efrat and Itai) 

Input : Graph G arising out of point sets A, B, and real number e > 0, and a 
matching M. 

Output : Layer graph C 


Li = exposed vertices of A ; 
i = l; 
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D = D t {B); 
repeat forever 

{ 

L 2f = 0; 

for each a € L 2 i-i 

while neighbor € (D, a) ^ 0 

{ 

b = neighbor* (£>, a); 
add b to L 2i ; 
delete* (D, b); 

} 

if L 2i is empty then { no augmenting path exists; Stop; } 
else 
{ 

if L 2i contains exposed vertices 

then { construction of C is complete; return £; } 

else Lx+i = all vertices of A adjacent to l 2i via edges of M; 

} 

i = 2 + l; 

} 


Next we describe the procedure for finding augmenting paths from the exposed 
vertices of Li to exposed vertices of L 2t . For this we need the following lemma. 

Lemma 4.1 (Efrat and Itai) Let M be a graph-matching of a bipartite graph G = 
(AUB, E), let TLbea set of edge-disjoint augmenting paths, and v be an intermediate 

vertex of some path of II. Then v cannot participate in any other augmenting path 
of n. 

To look for augmenting paths, D 2i s D t {Lu) is first constructed for each of the 
even layers L 2i C B of the layer graph C. Then a depth-first search is conducted, 
starting from an exposed vertex of L\. To advance from a vertex a € La- 1 > the 
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operation neighbor* (D 2 *, a) is invoked on the data structure £> 2 »• If this returns a 
vertex b € L 2 i then (a, b) is added to the current path and we advance to b. If 
neighbor* ( D 2 i, fl) returns 0, it indicates that a does not have any neighbors in L 2 i 
and hence does not lead to an exposed vertex of L 2t , so we should backtrack. To 
advance from b 6 La ( i < t ), if (6, a + ) € M then (fc, a + ) is added to the path 
and we advance from a + . Note that b is never exposed since all exposed vertices of 
B fi £ belong to L%t- If b € L^t is an exposed vertex, then an augmenting path is 
found. In this case M is increased and the vertices of this augmenting path from 
the appropriate LiS are deleted (because of Lemma 4.1), before starting to look for 
another augmenting path. To backtrack from a € £ 2 i-i (* > 2), if (b~, a) € M and 
a~ be the vertex preceeding b~ on the path found so far, then a and b~ are removed 
from the path and the search is continued from a~. If a 6 L x then it is simply 
deleted from L x . 

When no exposed vertices remain in L x , the search for augmenting paths in this 
particular phase is complete and we proceed to the next phase, where the layer graph 
C is again constructed and the whole process repeated. If during the construction of 
C one does not reach any exposed vertices of B, then it indicates that a maximum 
matching has been found. 

The time required to find all alternating paths in a single layer graph is 0(nT(n)). 
The time required to construct a layer graph was also shown to be 0(nT(n)). A 
theorem due to Hopcroft and Karp [HK73] states that Dinitz’s matching algorithm 
requires 0(%/n) phases. Hence Efrat and Itai’s scheme requires 0{n l - 5 T{n)) time. 

4.1.2 Faster Algorithms for Common Point Set Identifica- 
tion 

Making use of the scheme just described, in this subsection we show that replacing 
the Hopcroft and Karp’s algorithm in our approximation algorithms for common 
point set identification leads to an improvement in the overall running time. Recall 
that Algorithm 3.1 for deciding if there exists an a-LCP{A, B, e), had an indecision 
lateral of [! to («),lta(«)) “ d time complexity 0(n* s ) if .4 end B are point sets 
of cardinality O(n). We shall modify this algorithm to illustrate the improvement in 
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running time by using our new scheme. In that process it will also become clear that 
the same modification is extendable to all the other algorithms in a straightforward 
manner, to obtain the same speedup. 

First we describe the data structure D e (S) for our problem, that is made use of 
in Efrat and Itai’s algorithm. Arya et al. [AMN + 94] developed a data structure for 
answering nearest neighbor queries for a set of points in Given a set of n points 
S in and given any point q e 3? d , a nearest neighbor to q in S is any point 
s G S that minimizes d(s,q). For d > 3, there is no known algorithm for answering 
nearest neighbor queries that achieves both nearly linear space and polylogarithmic 
query time. The data structure proposed by Arya et al. reports in time O(logn), 
an approximated nearest neighbor of a query point q € 3R d . Approximate nearest 
neighbor is defined as a point s € S, such that for all s' € S , d(s, q ) < (1 + S)d(s f , q), 
where S > 0 is a predefined parameter. The construction of this data structure 
takes O(nlogn) time, O(n) space, and can also be dynamized so that deletion takes 
O(logn) time. 

Our new approximation algorithm for common point set identification uses the 
Efrat and Itai’s algorithm with Arya et aVs data structure for finding the maximum 
matching, instead of the Hopcroft and Karp’s algorithm as was done in the last 
chapter. Thus T(n), the upper bound on performing the operations neighbor, (S, q) 
and deleters', s) is O(logn). Recall from the last subsection that Efrat and Itai’s 
algorithm requires time 0(n u T(n)). Hence the new, but approximate matching 
algorithm has time complexity 0(n 15 log n) instead of 0(n 25 ). Let us refer to the 
new graph matching algorithm as AMQ (Approximate Matching), which takes as 
input two point sets A, B, and real numbers e > 0 and 5 > 0. To implement 
neighboring), Arya et aV s data structure is used to find a point s. If d(s,q) < 
(l + fy> then s is reported as neighbor, (S, q), otherwise neighbor, (5, q) is 0. We 
will show that the consequence of using this approximate matching in our decision 
algorithms is an increase in the indecision interval, along with a one sided error for 
a small range of values of the point location error c, both of which depend on S. 
This drawback is also reflected in the results obtained, by modifying this decision 
algorithm to approximate the size of a m ax (e), as we did in the last chapter. Recalling 
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the definitions of T P q from Definition 3.3, and G(T P q , e, A, B) from Definition 3.4, 
for triplets of points P and Q, we now state our algorithm and the consequent 
theorem. 

Algorithm 4.2 

Input : Point sets A, B, real numbers e > 0, 0 < a < 1, 5 > 0. 

Output : If a-LCP (A, B,e) exists then YES or DON’T KNOW, else NO or 
DON’T KNOW. 

for all triplets P = (Pii,Pj 2 ,Pi 3 ) from A 
for all triplets Q = (< In Ah Ah) f rom & 

{ 

compute T P q\ 

M = AM(T P q(A),B,€, <5); 

if ( M > o;min(|A|, |£|)) then return YES; 

} 

decision = NO; 

for all triplets P = {Pi^VinPh) f rorn A 
for all triplets Q - (q^q^q^) from B 

{ 

compute T P q ; 

M = AM{T pq {A),B, Se,sy, 

if (Af > Q!min(|A|, |£|)) then decision = YES; 

} 

if (decision == NO) then return NO else return DON'T KNOW; 

Theorem 4.2 Algorithm 4.2 always returns the correct answer about the existance 
of a-LCP(A,B,e) if 

(i) e> Seminia), or 

It either returns the correct answer or returns DON’T KNOW for values of e € 

it might return 
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any of the three possible answers - YES, NO, DON’T KNOW. A transformatio 
T pq , along with the bijective mapping l induced by the matching algorithm AM{\ 
that results in Algorithm 4-2 to return YES, correspond to a-LCP{A , B, (1 + <5)c). 

Proof: It follows from Lemma 3.1 that Algorithm 4.2 finds a transformation T PQ 

such that atleast omin(|A|, \B\) points of T PQ (A) are within 8e m i„(o:} distance of dis- 
tinct points of B. So the graph G(T PQ , e, A, B) always has a matching of size greater 
than or equal to amin(|A|, \B\) if e > 8e min (o;). However, the approximate matching 
algorithm AM ( ) considers a graph G 3 G(T P q, e, A, B) instead of G(T P q , e, A, B). 
This is because some edges of G are larger than e, but none are longer than (1 + 8)e 
as apparent from the discussion in the last subsection. AM ( ) thus finds a maximum 
matching in G where G(T P q, e, A, B) C G C G(T P q, (1 + 8)e, A , B). Any matching 
in G(T P q , e, A , B ) is also a matching in G, and if a matching in G(T P q, e, A , B) can 
be increased by an augmenting path then for a matching of the same size there 
also exists an augmenting path in G. Thus if G(T P q,c, A, B) has a matching of 
size k, then the approximate matching algorithm AM() finds a matching of size 
> k. Hence the algorithm returns YES for any e > 8e m j„(a). But since a maximum 
matching by AM{) might have some edges of length between € and (1 + <5}c, the 
labelling induced by AM{) along with T P q correspond to a-LCP{A , B, (1 -f 8)e). 

For any e < |e m t n (o:), no isometry enables in finding a-LCP(A,B, 8c). But 
since AM{ ) considers the graph G C G{T P q, 8(1 + <5)c, A, B), the indecision interval 
extends to on the left. If G(T P q, 8(1 + <5)e, A , B) does not contain a matching 
of size k, then G also does not have a matching of size k. For any e < there 

does not exist any isometry T PQ such that G(T PQ , 8(1 + 8)e , A, B) has a matching 
of size amin(|A|, |B|). Hence for such e Algorithm 4.2 always returns NO. 

For values of e e e m m (<*)), no isometry T P q enables in finding ct-LCP{A, 

B, e). But since AM{) considers a graph G C G(T PQ , (1 + 8)e, A, B), G might have 
a matching of size > amin(|A|, |B|). Hence the algorithm might return YES. The 
arguments for it returning NO or DON’T KNOW are exactly similar to those given 
for Theorem 3.2. 

For any e € [^i 7 ^ , there does not exist any transformation T P q for 

which the graph G(T P q,( 1 + 8)e,A,B) has a matching of size > amin(|A|, |B|). 
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Since AM() considers a graph G C G(Tpq, (1 + <5)e,.A, B), it can never return 
a matching of this size. Hence for such e the algorithm never returns YES. For 
e €E [€mm(oO> 8e m i n (o!)), it follows from Lemma 3.1 that the algorithm finds a trans- 
formation T pq , for which G(T PQ , 8e, A, B) has a matching of size > amin(|A|, |B|). 
For such a transformation AM{ ) also finds a matching of atleast this size. Hence 
the algorithm never returns NO for such values of e. | 

Theorem 4.3 If A and B are point sets of cardinality 0(n) then Algorithm J^.2 
runs in 0(n 7 - 5 log n) time. 

The corresponding algorithm that we gave in the last chapter, i.e. Algorithm 3.1, 
had a running time of 0(n 8 - 5 ). However, as spelled out by Theorem 4.2, the new algo- 
rithm suffers from the following shortcomings : (i) the indecision interval increases 
from [|€ m in(o:),8e m i n (a)) to U [e min (a),8e min (a)) (ii) although a- 

LCP(A, B, e) does not exist for any e € [— e m i n (a)) Algorithm 4.2 might return 
YES for such values of e. However, the corresponding Algorithm 3.1 never returns 
a wrong answer (iii) for any specified point location error e, if the algorithm returns 
YES then corresponding transformation along with the bijective mapping induced 
by the graph matching algorithm AM( ) correspond to a-LCP(A, B, (l+5)e). There 
was no such increase in point location error in the case of Algorithm 3.1. 

In Theorem 3.5 we saw that while trying to approximate Oi max (e), the lower bound 
on obtained by Algorithm 3.2 was < aw x (e) and on > max {a : e > 8e m i n (a)} (recall 
that if ai < <* 2 then e m ,- n (ai) < wM)* However, if we modify Algorithm 4.2 in 
the same way as was done for Algorithm 3.1, the lower bound on might be greater 
than owOe). But since the graph considered by AM() can have a matching of size 
atmost <w((l + <5)e) min(|A|, |B|), on will be < cw((l + 5)e). Hence we have 
max {a : e > 8c mi „(a)} < on < <w((l + <5)e). Using the same reasoning given in 
Theorem 4.2, a max (€) < < min {a : e < 

To reduce the indecision interval of [|e T ntn(o ; )>8 e min(Q:)) of Algorithm 3.1 to 
[|e min (a), 2e mi „(a)) in Algorithm 3.4, we used our exact algorithm for finding the 
largest common point set under pure rotation i.e. Algorithm 3.3. All the remain- 
ing algorithms of Chapter 3 hinge on this algorithm. Note that the graph match- 
ing algorithm is invoked within this exact algorithm. Replacing the Hopcroft and 
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Karp’s algorithm by the approximate matching algorithm AM() will reduce the 
running time of this algorithm from 0 (n 6 - 5 ) to 0(n 5 - 5 log n), and thereby speedup 
the overall running time of all the algorithms of the last chapter which make use 
of it. Recall from the last chapter that in this algorithm, given point sets A, B, a 
fixed point p, and a specified point location e, we constructed all possible domes 
D {j for pairs of points a, € A and bj € B. These domes partitioned the surface 
of the sphere S P around p, into a number of regions. Each point within a region 
corresponds to a rotation R such that if this region is formed by the intersection 
of the domes Aiju A a ja>'**> A*j*> then d(oj,, bj t ) < e for a*, € R(A), b Jt € B, 
l = 1,2 ,...,&. This gives rise to a bipartite graph where the nodes of the graph 
are points of A and B and the edges join pairs ( 04 , 64 ), l - 1,2, Clearly, 
every rotation within a region gives rise to the same graph. To apply our approx- 
imate matching algorithm AM() to find a maximum matching in such a graph, 
we require a rotation corresponding to any point in the region that gave rise to 
this graph. Note that in our sweep algorithm, that was used to traverse all the 
regions on the surface of S p , we constructed the bipartite graph corresponding to 
event points which were either intersection points, or the leftmost points of the 
domes. To replace the Hopcroft and Karp’s algorithm, we explicitly compute the 
rotation corresponding to each such event point, apply this rotation on the set A , 
and then use our approximate matching algorithm AM(). If R be the rotation 
found by Algorithm 3.3 that enables in finding a max {()-LCP{A,B ,€ ), under pure 
rotation, it implies that the graph G(R, e, A, B) has a maximum matching of size 
<Wc(e) min(|A|, |jB|), and it was found out by the Hopcroft and Karp’s algorithm. 
However, since our approximate graph matching algorithm AM ( ) considers a graph 
G, where G(R,e,A,B) C G C (?(J?, (1 + 5)e, A,B,), it might find a matching of 
size larger than the maximum matching in G{R,e,A,B). But for any rotation R, 
G(R, (1 + 6)e, A, B) can have a matching of size atmost + <5)e) min(jA|, |£j). 

So our new algorithm for finding the largest common point set under rotation re- 
turns a value a, where a maa .(e) < a < a mai ((l + 5)e). Use of this algorithm in 
Algorithm 3.4, has consequences similar to those seen in the case of Algorithm 4.2 
in the last two paragraphs. Recall that Algorithm 3.4 had an indecision interval of 
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^Cmtn (o;)) and time complexity of 0(n 8,5 ). Using our new algorithm for 
finding the largest common point set under rotation, changes the indecision inter- 
val to [ 2 {i+ 6 ) > "i+s ^ ^ [ e m<n(^)) 2e m ,- n (a)). Moreover the algorithm might err for 
values of e belonging to [-— 2 il_l ,e m i n (a)). However, this new algorithm runs in time 
0(n 7i5 logn). Using the same reasoning followed in the last two paragraphs, if this 
algorithm is made use of to approximate cw(e), the lower and upper bounds a, 
and a u will be as follows : max {a : e > 2e min (a)} < a, < a maa: ((l + 5)e) and 
^mai( £ ) <• &u min {a : e < Exactly similar results are obtained in the 

case of Algorithms 3.5 and 3.6 of the last chapter. 


4.2 Using Random Sampling 

In this section we show that application of standard random sampling techniques 
can reduce the time complexity of our algorithms by a considerable extent, how- 
ever, at the cost of a small failure probability. By now this technique has become 
fairly standard for this class of problems [ATT97, FKL+97, IR96]. We had briefly 
mentioned in Chapter 2 how this is done in the case of the alignment scheme, which 
is the basis of all our algorithms, except our exact algorithm for finding the largest 
common point set under rotation, i.e. Algorithm 3.3. The basic idea is that while 
searching for a transformation that enables in finding the largest common point set 
between two given point sets, instead of going through all possible transformations, 
only those arising out of a randomly sampled set of points are explored. If the size 
of the common point set is a considerable fraction of the size of the given point sets, 
then with a high probability we pick up points lying within the common point set, 
which thus enables in finding out a required transformation. Clearly, this technique 
can be directly applied to Algorithm 3.1, .where instead of going through all possi- 
ble triplets of points, only those corresponding to randomly sampled subsets of A 
and B are tested. However, recall from Theorem 3.2 that it was possible to bound 
the indecision interval of Algorithm 3.1, due to a transformation Tpq where P is 
a triplet that ‘enclosed’ the points of A belonging to the common point set. Since 
the algorithm went through all possible triplets of A, such a triplet was guaranteed 
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to be considered. We could not prove such a bound on the indecision interval in 
the randomized case. However for the remaining algorithms of Chapter 3, the ap- 
proximation ratios do hold in their corresponding randomized counterparts which 
we state below. We also bound the probability of failure in all the cases. 

Recall that all the algorithms of the last chapter except the first three, in the 
process of searching for a transformation that enables in finding the largest common 
point set, seperately computed the translation and the rotation components of the 
transformation. Towards this, for pairs of points a € A and b € B, if t a b denotes 
that translation that takes a to b, then all possible translations to* arising out of the 
point sets A and B were explored. For each such translation, our exact algorithm for 
finding out the largest common point set under rotation was invoked. Our random- 
ized algorithms are based on the scheme of exploring all translations corresponding 
to randomly sampled subsets of the given point sets, instead of the original ones. 
If the size of these randomly sampled subsets is smaller than that of the original 
sets, then this leads to fewer invocations of the algorithm for computing the rotation 
and hence reduces the overall running time. The speedup obtained depends on the 
ratio of the size of the original sets to that of the the sampled subsets. By carefully 
choosing the size of the subsets, we can bound the probability of failure of such 
a scheme. We illustrate this scheme by modifying Algorithm 3.4, and show how 
this affects the overall running time and failure probability. As we proceed, it will 
become clear that the same results will also hold in case of all the other algorithms. 

Our first use of randomization is in choosing a subset of the given point set A. 
Towards this we have the following algorithm and the consequent theorem. 

Algorithm 4.3 

Input : Point sets A, B, real numbers e > 0, Q < a < l. 

Output : Either YES or NO or DON’T KNOW. 

A' = a multiset of cardinality k, containing randomly sampled points of A; 
for each point a' € A! 
for each point be B 
{ 
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a' = Algorithm 3.3(t a > b (A),B,b,e)i 
if (a' > a) then return YES; 

} 

decision = NO; 
for each point a' € A' 
for each point b € B 

{ 

a' = Algorithm 3.3(t a >b(A), B,b,2e); 
if (a' > a) then decision = YES; 

} 

if (decision = = NO) then return NO else return DON’T KNOW; 


Before stating the theorem, let us assume that A and B are point sets of equal 
cardinality n. In the case of A and B being of unequal cardinality, the same type 
of analysis that we present here will follow. But this assumption leads to simplified 
expressions. 

Theorem 4.4 If point sets A and B are of cardinality n, and the cardinality of 
the randomly sampled multiset A' be a constant k, then Algorithm 4-3 runs in time 
0(n 7,i ). For any k > In j~], the algorithm returns YES with probability atleast 
q, for all c > 2e min (a). For € < §e m i„(o;) the algorithm always returns NO, and for 
2 c min(a) < € < € m »n(<*) it either returns NO or DON’T KNOW. 

Proof: Since Algorithm 3.3 has a time complexity of 0(n 6,5 ) it follows that 

Algorithm 4.3 runs in time 0(n 7-5 ). If \A' n a-LCP(A, B,e m i n (a))| > 1, then it 
follows from Theorem 3.8 that Algorithm 4.3 returns YES for all e > 2e m i„(a). We 
show that for Pr( | A'r\a-LCP{A, B,e min {a))\ > 1 ) > q to hold, it is sufficient that 
k > ^Inr*-. 

Ct 1 — 

Let a be a randomly sampled element of A. Then Pr( a G a-LCP(A, B, e m i n (a)) ) 
> a. Since A' is of cardinality k, Pr( A' n a-LCP{A, B, e min (o!) = 0) ) < (1 - <*)*. 
Hence for Pr( j A' n a-LCP{A, B, € mi „(a))| > 1 ) > <? to hold, it is sufficient that, 

1 - (1 - ot) k > q 
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Rearranging and taking logarithms of both sides, 

fcln(l— a) < ln(l— q) 

Expanding the logarithmic term, 

o? oft 1 

> h_ 

For this inequality to hold, it is sufficient that 

, . 1 , 1 

k > —In 

a 1 - q 

For any e < e m i n (a), no isometry can result in finding a-LCP(A t B,i). There- 
fore for such e the algorithm never returns YES. Moreover, for c < |e m m(«) no 
isometry enables in finding a-LCP(A, B, 2c). Thus such values of € always result in 
Algorithm 4.3 to return NO. | 

Our second randomized algorithm extends this concept of testing translations 
corresponding to only a randomly sampled subset of A. In this case we not only 
randomly sample a subset of A, but also randomly sample a subset of B and then 
test only the translations corresponding to these subsets. Rest of the algorithm is 
exactly similar to Algorithm 4.3. If A' and B' are randomly sampled multisets of A 
and B, then we invoke Algorithm 3.3 for all possible translations f«/y, d € A! and 
b' 6 B'. We thus have the following theorem. 

Theorem 4.5 If point sets A and B are of cardinality n, and the cardinality of the 
randomly sampled multisets A! and B' be constants k x and k 2 respectively, then the 
algorithm runs in time 0(n 6 5 ). For k x k 2 > any arbitrary constant c, the algorithm 
returns YES with probability atleast 1-e"^, for all e > 2e min (a). For e < |e m j n (a) 
the algorithm always returns NO, and for |e min (a) < € < € min (a) it either returns 
NO or DON’T KNOW. 

Proof: Let l be the bijective mapping underlying a-LCP(A,B,e min (a)). It 

follows from Theorem 3.8 that for any e > 2e min (a), the algorithm always returns 
YES if A 1 n a-LCP(A,B,e m i n (a)) ^ 0 and 1(a) € B', where a is any element of 
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A' n a-LCP(A, B,e min (a)). We show that for both the conditions to be satisfied 
with probability atloast 1 - e~ i n, it is sufficient that kik 2 > c. The rest of the proof 
is similar to that for Theorem 4.4. 

Let a be a randomly sampled element of A. Then Pr( a 6 a-LCP(A, B, e min (a)) ) 
> a. Assuming that a € a-LCP(A,B,e min (a)), let B x denote the event that the 
multiset B' contains 1(a). Then, Pr (Ei) = 1 - ( r ~) k2 . Let E 2 denote the event, that 
in the process of picking up a random element a € A and then randomly picking up 
k 2 elements of B with replacement, we never come across a pair a€ A,b e B such 
that o € a-LCP(A , B, c min (a)) and b = 1(a). Then Pr(B 2 ) < 1 - a {l - (^)* 2 }. 
Let £ 3 denote the event that in the process of randomly sampling ki elements of A 
with replacement, and for each such element randomly sampling k 2 elements of B 
with replacement, we come across atleast one pair a e A, b e B, such that a e a- 
LCP(A, B, € m in(a )) and 1(a) = b. This is precisely the event that we are looking for. 
Clearly, Pr(i? 3 ) >l-[l-a{l - We want this to be atleast 1 - e - ^ 

M^n] ~ e_ " 

This is equivalent to, 

-Prf}) * -? 


Expanding the In term gives 




For this to hold it is sufficient that 


-Mvfh? 


Expanding C 2 -*)* 2 gives 


k i 


ft 


h Mfo ~ ]) , \ > i 

2 n 2 J - n 


> 


coc 

n 


Hence, kik 2 > c. I 
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4.3 Practical Algorithms Exploiting the Stuctural 
Properties of Proteins 

All the algorithms presented so far have time complexity which arc relatively high 
degree polynomials of the size of the point sets. Recall from the last chapter that 
these 3-D point sets actually represent some drug or protein molecules, where each 
point is representative of an atom of the molecule. Since drug molecules are quite 
small in size, ranging from ten to a few hundred atoms, in practice these algorithms 
will be quite fast. However, since protein molecules are usually of the order of 
a few thousand atoms, in many practical situations such as database queries, the 
running time of our algorithms might be unacceptable. Till now we have not used 
the fact that our point sets are actually representative of some molecules. Since each 
atom of a molecule can be mapped only to a similar atom of the other molecule, this 
reduces the total number of possible transformations by a large extent. However, it is 
difficult to make use of this fact in the analysis of the running time of our algorithms 
in a general setting. In the case of particular molecules, by using a quantitative 
information about their structure, it will be possible to bound the exact running 
times apriori. In this section we show that by exploiting some general structural 
properties of proteins, it is possible to design algorithms with much smaller time 
complexities compared to those presented so far, however, at the cost of loosing the 
guaranteed performance bounds. But we expect that for all practical purpose they 
would perform quite well. Our algorithms are completely general and can be used 
for the common substructure identification in case of any protein molecules. 

4.3.1 Protein Structure 

A protein is composed of a chain of amino acids linked to each other by peptide 
bonds. There are three groups of atoms in each amino acid which constitute the 
backbone of the chain : the central atom C a , to which are attached on each side 
an N-H group and a carbonyl group C' = 0. The alkyl residue J ?, also bound to 
the C a , characterizes the nature of the amino acid but does not part in the 
backbone of the chain (see Figure 6). We will refer to the N-H group by the N 
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atom and the C' = O group by the C' atom. 


Peptide bond 



h o Rn+1 


Amino Acid n Amino Acid n+1 

Figure 6: Structure of a Protein Chain 


The sequence of amino acids fold in space to generate a complex three dimen- 
sional structure. Although there can be rotations around the C a -C' and C a -N 
bonds, and hence the geometry of the chain is weakly constrained, the geometry of 
the atoms attached to C a is perfectly determined. The three atoms N, C a and C 1 in 
each amino acid form a triangle which uniquely defines the position and orientation 
of the amino acid, and hence the entire protein structure, in space. Thus all the 
iV, C a and C' atoms act as a backbone or skeleton to which the alkyl residues R 
are attached. Since the C a -N and Ca-C bond lengths and the NC a C' bond angle 
are fixed, the skeletons corresponding to two common substructures of two proteins 
will be exactly congruent. Since correspondence between two triplets of points is 
sufficient to uniquely determine a rigid transformation in space, if we know the cor- 
respondence between two amino acids belonging to the common substructures of 
two protein molecules then we can compute the rigid transformation that results in 
the two skeletons to exactly coincide. We make use of this property in our algorithm 
in the next subsection. 
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4.3.2 An 0(n 45 ) Algorithm Using Protein Structure Infor- 
mation 

Given point sets A, B, and a real number c > 0, let l be the bijective mapping 
underlying a max {e)-LCP{A,B,e). Let p x € Of ma *(e)-LCP(A,£,<), pa be the point 
of oima x (t)-LCP{A,B,e) which is farthest from pi, and p 3 be the point of 
LCP(A,B,e ) such that the perpendicular distance of p 3 from the line pfpZ is maxi- 
mized. Now consider the isometric transformation I such that I{p\) = l(p\)\ I(pi), 
J(p 2 ) and l(p 2 ) are collinear; X{p x ), I(p 2 ), X(p 3 ) and l{p 3 ) are coplanar. Recall 
from Lemma 3.1 that the isometry I and the bijective mapping l correspond to 
a max (e)-LCP(A,B,8e). Note that the triplet (pi,p 2 ,p 3 ) defines a cylindrical region 
{P € 3ft 3 : |pip|. < \pm\ and dist {p,?ffi) < dist(p 3l (see Figure ?), where 
dist (p,pip 2 ) denotes the perpendicular distance of the point p from the straight line 
through pi and p 2 . 



Figure 7: The cylindrical region defined by the points pi,p 2 and p 3 

All points of amax(t)-LCP(A,B,e) lie on or within this cylindrical region and 
are atmost 8e distance away from their corresponding points in the set J 3, under 
the isometry I. In the case of our protein molecules, let Pi,p 2 ,p$ be either N, O a , 
or C' atoms lying on the skeleton or backbone of the common substructure and 
define the cylindrical region which enclose the entire skeleton of the substructure. 
Then the only parts of a max (e)-LCP(A,B,e) that can be outside this cylindrical 
region are some of the alkyl residues bound to the C a atoms lying on or within 



54 


Improvements Using Geometry, Randomization and Structure of the Molecules 


this cylindrical region. Hence the isometry I in this case takes each point of the 
common substructure to within 8e distance of its corresponding point in the set B, 
except possibly for some of the points corresponding to the alkyl residues which 
lie outside this cylindrical region. As mentioned in the last subsection, since the 
skeletons corresponding to the common substructures of two proteins are exactly 
congruent, any isometric transformation which maps three points on this skeleton 
to the corresponding points of the set B, will satisfy the properties of isometry I. 
Hence we have the following algorithm. 

Algorithm 4,4 

Input : Labelled point sets A and B corresponding to two protein molecules and a 
real number e > 0. 

Output : A real number a. 

M = 0; 

for each amino acid a € A 
for each amino acid b € B 

{ 

Tab = transformation that takes the N, C a and C' atoms of a to the corresponding 
atoms of b; 

M' — size of the maximum matching in G(T a b,Se,A,B)] 
if ( M ' > M ) then M = M'; 

} 

return Mf min (|A|, |BJ); 

Clearly, if the algorithm returns a, then there is a subset 5 of A of cardinality 
amin(|A|, jBj) which is 8f-congruent to some subset of B. If the skeleton or back- 
bone corresponding to the points of S, formed by the N , C a and C' atoms, enclose 
all the points of a m «,(e)-ICP(A, B, e), then it follows from the last two paragraphs 
that a > Qfma*(<0* However as already mentioned, since some of the alkyl residues 
of a max (e)~LCP(A, B, «) might be outside, it is possible that they will be more than 
8e distance away from their corresponding points in B under any of the transfor- 
mations tested by Algorithm 4.4. Under such circumstances a might be less than 
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a max (e). We expect that for most practical problems a > n maz (f) will hold, and 
even in cases where it does not, the difference will not be too large. Moreover, since 
any substructure identified by our scheme will ultimately be chemically cvcluated to 
test its effectiveness, the substructure found by our algorithm should only be large 
enough. If A, B are point sets of cardinality O(n), then transformations correspond- 
ing to 0(n 2 ) triplets are tested. The graph matching requires 0(n 3 5 ) time, Hence 
the overall complexity of the algorithm is 0(n 4,5 ). 

4.4 An 0(n 2 5 log n) Randomized Approximation Al- 
gorithm for Proteins 

In this section we give an algorithm, making use of the concepts presented in the 
last three sections. This algorithm is based on Algorithm 4 . 4 , with the approximate 
graph matching algorithm AM ( ) being used to find the maximum matching in the 
bipartite graph corresponding to each transformation. Secondly, instead of comput- 
ing transformations corresponding to all possible pairs of amino acids, we randomly 
sample a subset of amino acids from A and compute only those transformations 
corresponding to this subset, and the amino acids of B. 

Algorithm 4.5 

Input : Labelled point sets A and B corresponding to two protein molecules and 
real numbers e > 0, 5 > 0. 

Output : A real number a. 

A' — a randomly sampled multiset of k amino acids from A; 

M = 0; 

for each amino acid a' € A' 
for each amino acid b € B 
{ 

T 0 '& — transformation that takes the N 1 C a and C atoms of a 1 to the corresponding 
atoms of 6; 
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M' ~ AM(T a ib(A), B,8e t 5); 
if (A/' > M) then M = M'\ 

} 

return M/ min (|A|, |B|); 

Using the same reasoning as in Section 4.2, if the set A contains n amino acids, 
out of which atlcast m belong to cxmax(t)-LCP(A,B,t s) then the set A' contains 
atleast one of them with probability > q, for k > If this happens 

then Algorithm 4.4 computes atleast one transformation which maps the skeleton 
corresponding to the common substructure of A into the corresponding skeleton 
of B. Hence Algorithm 4.4 returns with probability atleast q, a real number a 
such that there exists a subset S C A of size amin(|A|, |B|) which is 8e(l + 5)- 
congruent to some subset of B. Moreover, in all practical problems we expect a to 
be greater than, or atleast very close to a maz (e). Since the set A' is of cardinality 
0(1), if there are n amino acids in the protein corresponding to the set B, then 
0(n) transformations are computed. Graph matching corresponding to each trans- 
formation takes 0(n l i logn) time. Hence the overall complexity of Algorithm 4.5 
is 0(n 2S logn). The corresponding algorithm of Akutsu returns a subset S C. A 
which is 8e-congruent to some subset of B and runs in time 0(n 8 ). However, it is 
guaranteed that \S\ > <Wc{*) nnn(jAl, jB|) will hold. 

It should be noted that instead of testing all the amino acids corresponding to 
the set B, we could have randomly sampled a subset of B in the same way as was 
done for set A, The resulting algorithm will have a time complexity of 0(n LS logn). 
However, the success probability would change from q to 1 - e~^ , where c is an 
arbitrary constant. 




Chapter 5 
Conclusions 


In this thesis we have presented several algorithms for identifying the structural 
similarities between two drug or protein molecules. An abstraction of this problem 
is that of finding the largest common point set under e-congruence, between two 
point sets in 3-D. This has important applications in biomolecular recognition and 
binding, synthetic drug design, and molecular database screening. Most commonly 
used algorithms which address this problem are based on character string comparison 
algorithms, and hence in the case of proteins consider only the primary structure 
of the protein chain. But all our algorithms directly address the inherent three 
dimensional structure of the molecules, which is essential to find out any m eaning ful 
similarity between two molecules. 

Although we have not presented any exact algorithm for the general problem, 
there is evidence that such an algorithm will have an exceedingly high running time. 
All our algorithms for the general isometry were approximation algorithms, with a 
guaranteed performance bound. Keeping in mind our primary application, given 
two molecules, since the allowable point location error e is fixed by the chemist, 
our objective has been to find out the largest common substructure between the 
two molecules which conform to this allowable error. In our abstract problem this 
amounts to finding out €«*«*(€). Towards this we have presented two classes of 
algorithms. The first approximates ^Ce) by satisfying the point location error 
e exactly, and the second approximately satisfies e. In both the cases we provide 
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guaranteed approximation ratios. 

Our first class of algorithms were based on the concept of approximate deci- 
sion algorithms for approximate congruence of point sets in 2-D, first introduced by 
Schirra [Sch92]. The first algorithm that we presented had an indecision interval of 
[|€ min (a), 8e mi „(a)). Then using our exact algorithm for finding the largest common 
point set under pure rotation, we reduced this indecision interval to [§e TO i n (o:), 2e m i n (a:) 
Both our algorithms had the same time complexity. However, from an implemen- 
tation point of view our first algorithm is considerably simpler. Although these 
indecision intervals might seem to be too large for these algorithms to be of practi- 
cal interest, it should be noted that they indicate the worst possible scenario. In all 
practical problems it is expected that we would obtain much better results. 

In the second class of algorithms, we have improved the approximation ratio 
obtained by Akutsu [Aku96], while keeping the time complexity same as that of his 
algorithm. Moreover, for the same approximation ratio, we have obtained algorithms 
with substantially better time complexities. 

We also presented two techniques, using which it was possible to reduce the 
running time of our algorithms. The first was an approximate graph matching 
algorithm, making use of the fact that the nodes of the graph were actually points 
in space, and the edges correspond to pairs of points which are within a specified 
distance of each other. Our second technique was the use of randomization. However 
both of these had some overhead. In the first case the indecision interval increased by 
a constant amount depending on the parameters of the graph matching algorithm. 
In the second case, as is common in all randomized monte carlo algorithms, our 
algorithm fails with some probability. Finally we showed one case, where making 
use of biological information can improve our algorithms by a large extent. 

We expect that for any problem, depending on the requirements it would be 
possible to put together techniques from those which we have presented, to obtain 
an algorithm which would work well in practice. 
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5.1 Future Research 

Throughout this thesis we have treated all molecules to be rigid bodies, and consid- 
ered only rigid transformations to superimpose one molecule on the other. Although 
such a treatment is adequate for comparing molecules with strong similarities, it will 
fail to identify weak similarities between pairs of molecules. So methods to overcome 
this need to be explored in future work. 

It is not sufficient to identify the similarities only between pairs of molecules. 
Most applications require the identification of the common substructure in a group 
of molecules. For m collections of n points on the real line, the largest common point 
set cannot be approximated to within an n e factor unless P = NP , and only weak 
positive results are known [AH]. Hence it is expected that extending our algorithms 
to find out the largest common point set of m point sets in 3-D, will be exceedingly 
difficult. Although Finn et al. [FKL + 97] have addressed this problem, no theoretical 
study was done. 

One can view our algorithms only as a basic paradigm. We have shown that 
using a very general structural property of proteins, it is possible to improve their 
time complexities by a considerable extent. In a similar manner, additional biolog- 
ical information can be incorporated into this basic framework. This emphasizes 
the importance of an interdisciplinary research involving Biology, Chemistry, and 
Computer Science, in this area. 
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